Problem 3 Consider model (7.6). Show that trajectories starting in the triangle T remain in T.Viral infections We first model the dynamics of a viral infection, such as hepatitis B or C, and are inter- ested in describing how the corresponding virus can spread and multiply in a person's body. Denote by X the average number of uninfected cells, which virions will try to in- fect; let Y be the average number of infected cells, and V be the average viral load (or the number of free virions in the body). Consider that uninfected cells are produced at a constant rate ) by the body, die at rate S X, and become infected at rate f(X, V) X, where f is some function of X and V. As a consequence, infected cells Y are created at rate f(X, V) X, and we assume they die at rate aY. Finally, free virions are pro- duced at a rate proportional to the number of infected cells KY, and are removed or destroyed at rate K V. If we consider that, as a first approximation, f is a linear func- tion of V, i.e. f(X, V) = bV, we have the following model.[!] dx = 1 - 8X- bVX, dt dy = b VX - aY, (7.1) dt dV =KY - KV, dt This model has six parameters and four variables. We can rescale time and V, but it is not a good idea to rescale X and Y independently since they both count cells and the term in bV X transfers cells from the X compartment to the Y compartment. We can therefore reduce Equations (7.1) to a model with three parameters. More precisely, leti=Cmvm, d7 d i = w my, (72) da- do 7 = _ 1, d'T y H 7 Akb where C = F, 17 = g, and p, = g are dimensionless parameters. We scaled time according to the death rate of normal cells. Alternatively, we could have scaled time according to the rate at which normal cells are produced by the body, i.e. we could have defined 1' = At/Xg, with X0 = 52/(kb). We could also have used a combination of these two time scales. In general, there is more than one possible way of defining dimensionless variables. The most convenient choice is often that which gives dimensionless parameters of order one, if at all possible. The classic endemic model For an endemic disease, births and deaths need to be taken into account, and the SIR model becomes ds dt = UN - us - as N ' dI I -ul + as- N - BI, (7.5) dt dR dt = -UR + BI, with initial conditions S(0) = So 2 0, I(0) = 10 > 0, and R(0) 2 0. Here the new parameters are the per capita death rate / and per capita birth rate v of the pop- ulation. By choosing / = v, the total population N = S + I + R is constant. In this case, using the same dimensionless variables as for the classic SIR model, we are left with a two-dimensional dynamical system, which reads, in dimensionless form, ds dT = n - ns - si, (7.6) di dT - (n+ 8)itsi, where n = v / a = M/a. As before, it is easy to check that trajectories starting in T remain in T (see exercises). The fixed points of (7.6) in the (s, i) plane areP1 = (1,0) and P2 || A + P" a A ._i | a | 5': V V The Jacobian of (E) is Jls'i) : (nii s (7+0) and JlP1)=(0n 1(nl+5))' Whether P2 is in T depends on the parameters 6 and 1]. More precisely, P2 6 'T 4:) 0 1, P1 is the only fixed point in 'T and since the eigenvalues of J(P1) are 1] and 1 (1] + 6), P1 is a stable node. All of the trajectories starting in T must converge to this xed point. which means that in the long run there are only susceptible individuals in the population. This is because those who have recovered from the disease eventually die and are replaced by new- borns, who are susceptible. This is illustrated in Figure 7.2, which shows the phase por- trait of (E), obtained with PPLANE, with 6 = 0.2 and 1] = 1. s =n-ns-si n = 1 i' = -(nto)i+si 6= 0.2 TTTTTTTTTTTTTTT TTTTTTTTTTTTTTTT .~ 0.5 0 0 0.5 S Figure 7.2. Phase plane of system (7.6), with 6 = 0.2 and n = 1, plotted with the software PPLANE. Only the dynamics inside T (not shaded) is relevant. If on the other hand n + 8 0, and R(0) 2 0. Here, S, I and R are the expected numbers of individuals in each compartment, and N = S + I + R is the total population. Note that births and deaths are not included in the model. This typically works for diseases which evolve over a short period of time, so that changes in the total population are negligible. In this model, I/N represents the fraction of infectious individuals. The product of this quantity with the contact rate a > 0 measures the average number of positive (i.e. giving rise to transmission of the disease) contacts per susceptible individual per unit of time. Since there are on average S susceptible individuals, the rate of change of S is aSI/N. The number of infected individuals increases by contact with susceptibles, and decreases due to recovery, at rate /SI, with B 2 0. By adding up the three equa- tions, one easily checks that dN/dt = 0, as expected since we neglected births and deaths. System (E) can thus be reduced to a two-dimensional system of ordinary dif- ferential equations, by omitting the last equation for R. The remaining two equations may be written in dimensionless form by letting S _ I s=i,z=i,and7'=at. N N Then, the first two equations of (Q) become ds _ 7 2 .5'1 dT ' (7.4) dz' 7 = 31' 6i d'T ' where 6 = /a > 0. The quantity :7 = 1/15 is the contact rate a multiplied by the characteristic time 1/ during which a person remains infectious. It is called the con- tact number of the disease, and in this case is equal to the basic reproduction number R0 of the infection described by the SIR model. Since (1 s i)N = R Z 0, Equations (M) only make biological sense ifs and 1' remain positive and such that s + 'i S 1, provided initial conditions satisfy these re- quirements. In other words, trajectories of the dynamical system (M) that start in the triangle T={(s,i)|520,i20,s+i1} should remain in T. To check this, consider the dynamics on the boundary of T. First assume that s = 0 and i S 1. Then dS/d'T : 0, i.e. 5 remains equal to zero, and 1' de- creases towards zero, but does not become negative. Similarly, ifi : 0, then both 3 and i remain constant (what is the biological significance of this fact?). Finally, if d s +z' = 1. then (17\" +1") = 6i 3 0, so that s +z' will not increase past the '7' value 1. System (Q) has an infinite number of fixed points in T, which are such that i : 0 and s is arhitrarv , P OneDrive Screenshot saved Figure Z1. Phase plane ofsystem {Z_4), with 6 = 0.2, plotted with the software PPLANE. Only the dynamics inside T (not shaded? is relevant Figure 7.1 shows the phase portrait of (H), obtained with PPLANE. with 6 = 0.2. In this case, we see that all trajectories in T converge towards one of the fixed points, i.e. limi = 0. This means that the epidemic eventually dies out, and we are only left t>oo with susceptible individuals and/or those who have recovered from the disease. Summary This chapter illustrates how the modeling concepts discussed in the previous chapters may be applied to the description, in terms of ordinary differential equations, of the dy- namics of infectious diseases and the spread of epidemics. Moreover, it provides the reader with a basic introduction to the terminology and applications of compartmental models. It should be clear by now that models of arbitrary complexity may be built from the simple tools discussed in this text. The modeling process is always the same, no matter how involved the model is. The methods of analysis, in terms of maps or differential equations, are also similar, but become more complicated as the dimension of the model is increased. In particular, three-dimensional continuously differentiable dynami- cal systems may exhibit chaos, the understanding of which requires more advanced techniques than those discussed here. The section on further reading includes texts on dynamical systems and chaos that the reader may want to consult