Fricative.math.ttu.edu

PROJECT 1: Epidemic Model: To Treat or Not to Treat? Topics: Nonlinear differential equations, equilibria, stability, basic reproduction number, continuous-time Markovchain model, drug resistance.
Antiviral drug resistance is an important concern in public health. The emergence and spread of drugresistance leads to treatment failure, limits the effectiveness of subsequent treatments, and resultsin the evolution of more virulent strains. Antiviral drugs such as amantadine and rimantadine (M2protein inhibitors) have been used for influenza A prevention and control. Such treatments are about80% to 90% effective at preventing infection, and can reduce the duration of the flu by about one tofive days, if given within 48 hours of infection. However, drug-treated persons infected with influenzamay shed resistant viruses within five to seven days after treatment. Recently, high levels of resistanceto amantadine and rimantadine have been detected among influenza A viruses. It has been shownthat resistance to amantadine and rimantadine increased from 0.4% in 1994-1995 to 12.3% in 2003-2004 among circulating influenza viruses collected worldwide between October 1994 and March 2005[4]. Neither of these drugs are recommended for swine influenza A (H1N1) prevention [4, 5]. Otherantiviral drugs, such as the neuraminidase inhibitors zanamivir (trade name Relenza R , inhalant) andoseltamivir (trade name Tamiflu R , oral medication), are currently recommended [4]. The viral enzymeneuraminidase functions to detach newly-produced virus from the host cell while the neuraminidaseinhibitors prevent the release of newly-produced virus from host cell membranes. The drugs are mosteffective when given early in the disease process when viral replication is high[6]. Both zanamivir andoseltamivir are recommended for antiviral chemoprophylaxis.
To model drug resistance during an epidemic, we first model the epidemic process using a standard SIR epidemic model. The population is divided into individuals that are susceptible (healthy), infectious, and recovered, denotedas S, I, and R, respectively. A compartmental diagram of the SIR model is given below. The transmission parameteris β and recovery parameter is γ (1/γ is the average length of infectivity). The total population size is N . Whether an outbreak occurs depends on the basic reproduction number, R0 = The differential equations that describe the dynamics of the SIR epidemic model are This simple model does not include a latent period, the time after exposure until a person begins shedding thevirus. Replication occurs in the epithelial tissue of the lungs. To formulate a continuous-time Markov chain modeleach of the events and their corresponding probabilities must be identified. In particular, let the changes in theprocess be denoted as (∆S, ∆I) = (S(t + ∆t) − S(t), ∆I(t + ∆t) − I(t)). Then the following table gives the eventsand their associated probabilities. Notice that we did not include R(t) because none of the events depend on R(t).
In a Markov chain, there is an exponential amount of time between each event. The Markov property meansthat each event only depends on the present, not the past history, because the exponential distribution has thememoryless property [1]. The exponential distribution for the time between events given the size of the susceptible and infective populations are S(t) and I(t) has a parameter λ = β function for the exponential distribution for the interevent time T is F (t) = 1 − eλt = Prob{T ≤ t}. A uniformrandom variable U can be used to compute the interevent time T [1], Prob{F −1(U ) ≤ t} = Prob{F (F −1(U )) ≤ F (t)} = Prob{U ≤ F (t)} = F (t), A Matlab program that numerically computes the solution to the ODE model and generates sample paths of the continuous-time Markov chain (CTMC) model is given at the end of this project (STOCHSIRProject.m).
Generating many sample paths provides an estimate of the probability there is an epidemic or an outbreak, i.e.,whether the number of cases is sufficiently large (> 20). It can be shown for the SIR model that an approximateprobability of an outbreak is provided R0 > 1 and where k is the initial number of infected individuals [1]. The probability of an outbreakincreases as R0 or k increase. This relationship holds only for the SIR, SIS or simple stochastic epidemic modelsand not for more complex models [1].
Exercise Use the Matlab program STOCHSIRProject.m and verify that the probability of an outbreak (fork = 1, . . . , 10 and R0 = 2, 3, 4) agrees with the approximation (1).
It is also possible to calculate the distribution for the time until the epidemic ends. This can be done in thesimulation by keeping track of the first time t such that I(t) = 0. This state is known as an absorbing state inMarkov chain theory because once I(t) = 0 it stays at zero.
How do we add drug treatment? How does the model change? It depends on when treatment is given.
Treatment can be given to susceptible individuals as a preventive measure known as chemoprophylaxis or given toinfected individuals to alleviate their symptoms. Here, we only consider chemoprophylaxis. Antiviral drugs maydecrease transmission, theseverity of symptoms, or the length of infectivity.
mutate, mutant virus strains that are able to escape the action of the drug, resistant strains, may become dominant.
The mutant virus strain may no longer be sensitive to the drug treatment. Thus, we consider a new model for thepopulation that is under drug treatment. Let Spr =susceptible individuals taking drugs prophylactically, Itr =infected treated individuals, and Ir,tr = infected treated individuals who shed drug-resistant virus. This is asimplified model because it does not include a latent period and the possibility of multiple resistant strains. Thenew compartmental diagram is given below.
The differential equations corresponding to a new model with drug resistance are as follows: New parameters are βr, r, κ, and p. Transmission may be reduced if infected by a resistant strain, βr ≤ β or ifan individual has been treated, pβ ≤ β. In addition, treated individuals who have become infected may have ashorter period of infectivity, 1/(rγ) = f /γ, f = 1/r < 1. The most important parameter is κ, resistance developsin treated individuals at a rate κ. For the continuous-time Markov chain (CTMC) model, there are five differentevents affecting Spr, Itr and Ir,tr. The changes and the corresponding probabilities are given in the table below.
Although this ODE model with treatment and drug resistance is more complex than the original SIR model, it is easy to check that the total population size remains constant. The sum of all differential equations is zero,so that the total population size is N . We will compute a treatment reproduction number, RT , using the nextgeneration matrix approach. The disease-free equilibrium (DFE) is ¯ Exercise Check by substituting the DFE into the differential equations, that the DFE is a constant solution.
To calculate the treatment reproduction number, only consider the infectives Itr and Ir,tr. We will call it treatment reproduction number to distinguish it from the basic reproduction number for the SIR model. The ratesof new infections are given in the vector F and the rates of transitions out of the two infective groups are given inthe vector V, where Exercise: Calculate the Jacobian matrix F and V for each of the vectors F and V with respect to x = Itr andy = Ir,tr and evaluate at the DFE. Then calculate F V −1 and show that the spectral radius RT = ρ(F V −1) is and evaluate at the DFE. Find the eigenvalues of F V −1 (they are both positive). The largest one is the treatmentreproductive number. There are two ways the disease is spread; this is why RT is the maximum of two values.
The first ratio is from the treated individuals and the second ratio is from the resistant individuals. Whicheverratio is greatest determines whether an outbreak occurs, whether RT > 1. Is it possible for individuals infectedwith the resistant strain to cause outbreaks. Compare RT with R0. Can RT be greater than R0? Some basic parameter values are given in the following table for the Matlab program STOCHDrugResist- Proj.m The population size is N = 1000. The recovery rate is γ = 0.25 which means the average length of theinfective period is 1/γ = 4 days (without treatment or for a resistant strain). The parameter f = 1/r shortens theinfectious period for the treated individuals by a factor ranging from 0.5 to 1, f /γ. If f = 0.5, then the infectiousperiod for the treated individuals is f /γ = 2 days. The transmission rate of the treated individuals can be reducedby a factor p, where p ranges from 0.1 to 1.0. In addition, the initial number of infected and treated individualswill be important in determining whether an outbreak occurs, Itr(0). The initial number ranges from 1 to 10individuals.
Exercise: Select different values of R0, p, f , and Itr(0) and compute the probability that resistance develops,the probability of an outbreak, etc. Discuss your results. Change the Matlab programs and make some otherassumptions about the dynamics of drug resistance. What happens if βr = β/2 or if βr = 2β, that is, the resistantstrain is less or more transmissible than the sensitive strain? Check the literature on drug resistance for someother ideas (some references below [2, 3, 7]).
Table 4: Probability that Drug Resistance Develops.
%Deterministic and Stochastic SIR Epidemic Models%Length of infectivity 1/g=4 days%Population size is 1000. Time is 120 days or 4 monthsclear allg=.25; N=1000; % Parameter Valuesdt=.05; tim=120; time=tim/dt;sim=1000;chk1=1;chk2=1;while chk1==1 R0=input(’ Enter R0, Between 0.5 and 5.0: ’); %Choose R0disp(’ ’);chk1=2; ii=input(’ Enter Initial I, Between 1 and 10: ’); %Choose Initial Idisp(’ ’);chk2=2; if( ii<1 | ii>10 ) chk2=1; endendii=round(ii);i(1)=ii; s(1)=N-ii; counti=ii;for tt=1:time % Euler’s method i(tt+1)=i(tt)+dt*((beta/N)*i(tt)*s(tt)-g*i(tt));s(tt+1)=s(tt)+dt*(-(beta/N)*i(tt)*s(tt)); endTotalCases_ODE=round(counti) % Total Cases after 120 daysset(0,’DefaultAxesFontSize’,16)set(gca,’fontsize’,16);plot([0:dt:tim],i,’k--’,’linewidth’,2);hold onEpid=0;for k=1:sim % 1000 sample paths t(1)=0; i(1)=ii; s(1)=N-ii;j=1;while i(j)>0.5 & t(j)<time ev1=beta*i(j)*s(j)/N;ev2=g*i(j);prob1=ev1;tot=ev1+ev2;t(j+1)=-log(u1)/tot+t(j); if max(i)>20; Epid=Epid+1; endif k==1 stairs(t,i,’r-’,’linewidth’,2); endif k==2 stairs(t,i,’g-’,’linewidth’,2); endif k==3 stairs(t,i,’b-’,’linewidth’,2); endhold on endhold offxlabel(’Days’); ylabel(’Number of Cases’);legend(’ODE’,’Stoch_1’,’Stoch_2’,’Stoch_3’)ProbOutbreak_CTMC=Epid/sim %Probability Outbreak (> 20 cases) %Deterministic and Stochastic Models with Drug Resistance%Length of infectivity 1/g=4 days%Rate of conversion of nonresistant to resistant k=.005%Population size is 1000. Time is 120 days or 4 monthsclear allN=1000; g=0.25; k=0.005;time=120; dt=0.05; tim=time/dt; sim=1000;chk0=1;chk1=1;chk2=1; chk3=1;while chk0==1 R0=input(’Enter R0, Between 0.5 and 5: ’); % Choose R0disp(’ ’); chk0=2; if( R0<0.5 | R0>5. ) chk0=1; endendbeta=R0*(g);betar=beta; % Transmission of Resistant Strains=NonResistant%Choose the effectiveness of treatment on%reduction of transmission: pwhile chk1==1 p=input(’Enter p, Between 0.05 and 1: ’);disp(’ ’);chk1=2; if( p<.05 | p>1. ) chk1=1; endend%Choose effectiveness of treatment on%reduction of length of infectivity by a factor f: f/gwhile chk2==1 f=input(’Enter f, Between 0.5 and 1: ’);disp(’ ’);chk2=2; if( f<0.5| f>1) chk2=1; endendr=1/f;%Calculation of R0 and RTR0=beta/gRT=max(p*beta/(r*g+k),betar/g)%Choose initial number of treated that are infected.
%Number of resistant cases are set to zero.
while chk3==1 ii=input(’Enter Initial Itr, Between 1 and 10: ’);disp(’ ’);chk3=2; if( ii<1 | ii>10 ) chk3=1; endendii=round(ii);Spr(1)=N-ii; Itr(1)=ii; Irtr(1)=0;countitr=ii; countirtr=0;for t=1:tim f1=-(p*beta*Itr(t)+betar*Irtr(t))*Spr(t)/N;f2=(p*beta*Itr(t))*Spr(t)/N-(r*g+k)*Itr(t);f3=(betar*Irtr(t))*Spr(t)/N-g*Irtr(t)+k*Itr(t);countitr=countitr+dt*(p*beta*Itr(t)*Spr(t)/N); countirtr=countirtr+dt*(betar*Irtr(t)*Spr(t)/N+k*Itr(t)); Spr(t+1)=Spr(t)+dt*f1;Itr(t+1)=Itr(t)+dt*f2;Irtr(t+1)=Irtr(t)+dt*f3; endTotalResist_ODE=round(countirtr) %Total Resistant Cases (ODE)TotalNonResist_ODE=round(countitr) % Total NonResistant Cases (ODE)set(0,’DefaultAxesFontSize’, 16); plot([0:dt:time],Itr,’g-’,[0:dt:time],Irtr,’r-’,’linewidth’,2);legend(’I_{tr}’,’I_{r,tr}’);xlabel(’Days’); ylabel(’Number of Cases’);hold onscountitr=ii; scountirtr=0.;j=1;spr(1)=N-ii; itr(1)=ii; irtr(1)=0; tt(1)=0; while tt(j)<time & itr(j)+irtr(j)>0.5; % One Sample Path u1=rand;u2=rand;ev1=p*beta*itr(j)*spr(j)/N;ev2=betar*irtr(j)*spr(j)/N;ev3=r*g*itr(j);ev4=g*irtr(j);ev5=k*itr(j);prob1=ev1;prob2=prob1+ev2;prob3=prob2+ev3;prob4=prob3+ev4;tot=prob4+ev5;tt(j+1)=tt(j)-log(u1)/tot;if u2<=prob1/tot;spr(j+1)=spr(j)-1;itr(j+1)=itr(j)+1;irtr(j+1)=irtr(j); scountitr=scountitr+1;elseif u2>prob1/tot & u2<=prob2/tot;spr(j+1)=spr(j)-1;itr(j+1)=itr(j);irtr(j+1)=irtr(j)+1;scountirtr=scountirtr+1;elseif u2> prob2/tot & u2<=prob3/tot;spr(j+1)=spr(j);itr(j+1)=itr(j)-1;irtr(j+1)=irtr(j);elseif u2>prob3/tot & u2<=prob4/tot;spr(j+1)=spr(j);itr(j+1)=itr(j);irtr(j+1)=irtr(j)-1;elsespr(j+1)=spr(j);itr(j+1)=itr(j)-1;irtr(j+1)=irtr(j)+1;scountirtr=scountirtr+1;endj=j+1; endplot(tt,itr,’g--’,tt,irtr,’r--’,’linewidth’,2);hold offTotalResist_CTMC=scountirtr %Total Resistant Cases (One Sample Path)TotalNonResist_CTMC=scountitr %Total NonResistant Cases (One Sample Path)Resist=0; j=1;scountitr=ii;scountirtr=0;spr(1)=N-ii; itr(1)=ii; irtr(1)=0; tt(1)=0; while tt(j)<time & itr(j)+irtr(j)>0.5; u1=rand;u2=rand;ev1=p*beta*itr(j)*spr(j)/N;ev2=betar*irtr(j)*spr(j)/N;ev3=r*g*itr(j);ev4=g*irtr(j);ev5=k*itr(j);prob1=ev1;prob2=prob1+ev2;prob3=prob2+ev3;prob4=prob3+ev4;tot=prob4+ev5;tt(j+1)=tt(j)-log(u1)/tot;if u2<=prob1/tot;spr(j+1)=spr(j)-1;itr(j+1)=itr(j)+1;irtr(j+1)=irtr(j);scountitr=scountitr+1;elseif u2>prob1/tot & u2<=prob2/tot;spr(j+1)=spr(j)-1;itr(j+1)=itr(j);irtr(j+1)=irtr(j)+1;scountirtr=scountirtr+1;elseif u2> prob2/tot & u2<=prob3/tot;spr(j+1)=spr(j);itr(j+1)=itr(j)-1;irtr(j+1)=irtr(j);elseif u2>prob3/tot & u2<=prob4/tot; spr(j+1)=spr(j);itr(j+1)=itr(j);irtr(j+1)=irtr(j)-1;elsespr(j+1)=spr(j);itr(j+1)=itr(j)-1;irtr(j+1)=irtr(j)+1;scountirtr=scountirtr+1;endj=j+1; endif scountirtr>10 Resist=Resist+1; end endProbResistOutbreak=Resist/sim % Probability Resistant cases >10 [1] Allen, L. J. S. 2003. An Introduction to Stochastic Processes with Applications to Biology. Prentice Hall, Upper [2] Blower, S., E. Bodine, J. Kahn, and W. McFarland. 2005. The antiretroviral rollout and drug-resistant HIV in Africa: insights from empirical data and the theoretical models. AIDS 19: 1-14.
[3] Blower, S. M. and T. Chou. 2004. Modeling the emergence of the ‘hot zones’: tuberculosis and the amplification dynamics of drug resistance. Nature Medicine 10: 1111-1116.
[4] CDC. 2006. High levels of amantadine resistance among influenza H3N2 viruses and interim guidelines for use of antiviral agents - United States, 2005-06 influenza season. MMWR Jan. 17, 2006, 55: 44-46.
[5] CDC. 2009. Antiviral drugs and H1N1 flu (swine flu). http://www.cdc.gov/H1N1flu/antiviral.htm [Accessed: [6] Bell, E. A. 2009. Pharmacotherapy for novel influenza A (H1N1). Infectious Disease News June 1, 2009.
http://www.infectiousdiseasenews.com/article/40914.aspx [Accessed: June 30, 2009] [7] Rong, L., H. Dahari, R.M. Ribeiro, and A. S. Perelson. 2010. Rapid emergence of protease inhibitor resistance in Hepatitis C virus. Science Translational Medicine [online May 26,2010].
[8] Stilianakis, N. I., A. S. Perelson, and F. G. Hayden. 1998. Emergence of drug resistance during an influenza epidemic: insights from a mathematical model. J. Infect. Dis. 177: 863-873.

Source: http://fricative.math.ttu.edu/past/MAAPREP/2010/pix-papers/Project1.pdf