BIOS13:SI dynamics

SI dynamics (8p)

Consider the spread of an incurable disease in a population, i.e. a disease from which there is no escape. Its spread can be modeled with an SI model, denoting S the number of susceptible and I the number of infected. Assume the dynamics of S and I follow:
$$
\begin{cases}

\cfrac{dS}{dt}=bS(1-\theta S)-\beta SI-\mu_S S\\

\cfrac{dI}{dt}=\beta SI -\mu_I

I
\end{cases}
$$
where b is the per capita birth rate (at low density), θ is the strength of density dependence, controlling the birth rate, 𝛽 is the ‘contact rate’, 𝜇_S is the death rate
(mortality) of the susceptible and 𝜇_I is the death rate of infected. It is thus assumed
that infected individuals are too sick to reproduce (or even affect the density dependence). All model parameters are positive constants.

Initially, ignore the disease (set I = 0) and consider the dynamics of the susceptible
population on its own.

• Without the disease, what is the equilibrium population size (S )? (1p)*
  $$f(S)=\frac{dS}{dt}=bS(1-\theta S)-{\mu }_{S}S=S(b-b\theta S-{\mu }_S)=0$$

  $$S_1=0S_2=\frac{b-{\mu }_S}{b\theta }$$

$$S^{\ast }=\frac{b-{\mu }_S}{b\theta }$$

• For what parameter values is S > 0? (1p)*

  Because all model parameters are positive constants, when b > 𝜇_S, S* > 0

• Show that this S is stable. (1p)*
  $$f^{\prime }(S^{\ast })=b-2b\theta S-{\mu }_S={\mu }_S-b$$
  We have already assumed that b > 𝜇_S, so $f^{\prime }(S^{\ast })<0$, the S* is stable.

Now let’s consider the full dynamics, i.e. let I > 0.

• Write a script in R that simulates the SI model and plots the resulting S and I over time. Use parameters 𝛽 = 0.01, 𝑏 = 1, 𝜃 = 0.001, 𝜇_S = 0.5, 𝜇_I = 2 and run the simulation 100 time units, starting at 𝑆 = 400, 𝐼 = 1. (1p)

  # plots the resulting S and I over time
  SI_sys = function(t, SI, p){
    S = SI[1]
    I = SI[2]
    dsdt = p$b*S*(1-p$theta*S)-p$beta*S*I-p$mu_s*S
    didt = p$beta*S*I - p$mu_i*I
    list(c(dsdt, didt))
  }
  
  timeVec = seq(1, 100, by=0.1)
  p = list(beta = 0.01, b = 1, theta = 0.001, mu_s = 0.5, mu_i = 2)
  SI = c(S = 400, I = 1)
  out = ode(y = SI, func = SI_sys, times = timeVec, parms = p)
  time = out[,'time']
  S = out[,'S']
  I = out[,'I']
  
  plot(time, S, type = 'l', xlab = 'Time', ylab = 'S&I', col = 'blue', ylim = c(0, max(S)))
  lines(time, I, type = 'l', col = 'red')
  legend("topright",  c("S","I"), lty=c(1, 1), col=c("blue", "red"))
  

  [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-eJbO7HbA-

• Extend the script to also plot the simulation in the SI phase plane. (1p)

  See below.

• Extend the script again to plot the isoclines of S and I in the phase plane. (2p)

$$\frac{dI}{dt}=0\Rightarrow S^{\ast }=\frac{{\mu }_I}{\beta }=200$$

$$\frac{dS}{dt}=0\Rightarrow I=\frac{b-b\theta S-{\mu }_S}{\beta }$$

Given S^*=200: I^*=30

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-xWkIk1MV-

  SI_isoclines = function(p){
    se = p$mu_i/p$beta
    s = c(0,500)
    i = (p$b-p$b*p$theta*s-p$mu_s)/p$beta
    plot(s,i, type = 'l', col = 'red', ylim = c(0, 150), xlab = 'S', ylab = 'I')
    lines(c(se, se), c(0, 150), col = 'blue')
    legend("topright",  c("S isocline","I isocline"), lty=c(1, 1), col=c("blue", "red"))
  }
  SI_isoclines(p)
  lines(S, I)
  points(S[1], I[1], pch = 0)

• What is the new equilibrium population size (S+I)? You can solve this either by simulations or calculations. However you solve it, write an R script that plots the equilibrium total population as a function of 𝜇_I, where 𝜇_I ranges from 1.0 to 4.0. The other parameters should take values as above. (1p)
  $$f({\mu }_I)=S^{\ast }+I^{\ast }=\frac{{\mu }_I-b-b\theta \frac{{\mu }_I}{\beta }-{\mu }_S}{\beta }=90{\mu }_I+50$$

  mu_i2 = c(1.0, 4.0)
  N = mu_i2/p$beta+(p$b-p$b*p$theta*(mu_i2/p$beta)-p$mu_s)/p$beta
  plot(mu_i2, N, type = 'l', xlab = 'mu', ylab = 'Total population')
  

  [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-F5vwIopi-