 24-09-2026            >eR-BioStat and STAEM-3D
Modeling infectious diseases using R - a multiple populations model
Ziv Shkedy and Rashider Aloni based on Chapter 2-4 in the book Modeling infectious diseases parameters based on serological and social contact data (https://www.simid.be/modeling-infectious-disease-parameters-based-on-serological-and-social-contact-data/)
## load/install libraries
.libPaths(c("./Rpackages",.libPaths()))
library(knitr)
library(tidyverse)
library(deSolve)
library(minpack.lm)
library(ggpubr)
library(readxl)
library(gamlss)
library(data.table)
library(grid)
library(png)
library(nlme)
library(gridExtra)
library(mvtnorm)
library(e1071)
library(lattice)
library(ggplot2)
library(dslabs)
library(NHANES)
library(plyr)
library(dplyr)
library(nasaweather)
library(ggplot2)
library(gganimate)
library(av)
library(gifski)
library(foreach)
library("DAAG")
library(DT)
library(TeachingDemos)
library(gridExtra)1. The SIR model with multiple populations
Low of mass action in one population setting
The law of mass action discussed before states that the number of newly infected individuals at time \(t\) depends on the number of infected individuals (\(I\)), the number of susceptibles (\(S\)) and the transmission rate between these two groups (\(\beta\)). In its simplest form, the mass-action principle states: \[ \mbox{Number of new cases} = \beta I \times S = \lambda S. \] The underlying assumption behind the latter equation is that infected and susceptible individuals mix in a random manner (homogeneous mixing) and therefore \(\beta\) is age and time independent.
Multiple populations
In this tutorial we discuss transmission settings in which the social structures in the population determine the mixing patterns in the population. We assume that the population is constructed from multiple sub-populations which may or may not interact with each other. As a result, the transmission process of the infection is governed by the mixing patterns in the population.
Heterogenous mixing patterns and Sub-populations
The above assumption about random mixing in the population usually does not hold. Most populations do not mix in a random fashion but are constructed from sub-populations in and between which individuals mix such as for example: different age groups within a school, contacts within households and sexual contacts within the population. Capasso (2008) refers to this mixing pattern in the population as a positive feedback epidemic system. A central characteristic of a positive feedback epidemic system is the mixing or Who Acquires Infection From Whom (WAIFW) matrix which governs the transmission process. For example, considering a population with two sub-populations A and B, the mixing pattern between these two groups can be describe by the following WAIFW matrix: \[ C= \left( \begin{array}{cc} \beta_{aa} & \beta_{ab}\\ \beta_{ba} & \beta_{bb} \end{array} \right). \] The mixing rates for individuals from group A with those of group B is denoted as \(\beta_{ab}\). Individuals can also mix with individuals from the same group, \(\beta_{aa}\) and \(\beta_{bb}\) denote the within-group mixing rates for groups A and B, respectively.
2. Example 1: a SIS model for Gonorrhea
Transmission model for Gonorrhea
Gonorrhea is a disease transmitted through the population by heterosexual contacts. Since immunity to reinfection does not exist, infected individuals, after recovering, move back to the susceptible class. Capasso (2008) proposed a SIS transmission model for Gonorrhea in which the two subpopulations, males and females interact according to the following mixing matrix
\[ C= \left( \begin{array}{cc} 0 & \beta_{fm}\\ \beta_{mf} & 0 \end{array} \right). \]
The SIS transmission model for gonorrhea can then be described by the following system of (ordinary) differential equations: \[ \left \{ \begin{array}{l} \frac{dS_f}{dt} = - \beta_{fm} S_f(t) I_m(t)+ v_f I_f(t),\\[2ex] \frac{dI_f}{dt} = \beta_{fm} S_f(t) I_m(t)- v_f I_f(t),\\[2ex] \frac{dS_m}{dt} = - \beta_{mf} S_m(t) I_f(t)+ v_m I_m(t),\\[2ex] \frac{dI_m}{dt} = \beta_{mf} S_m(t) I_f(t)- v_m I_m(t). \end{array} \right. \]
Here, \(v_f\) and \(v_m\) are the recovery rates for females and males, respectively. Note that the transmission model assumes that mixing between sub-populations occurs in a randomly fashion and therefore the number of new cases equals \(\beta_{fm} S_f(t) I_m(t)\) and \(\beta_{mf} S_m(t) I_f(t)\) for females and males, respectively. We further assume that \(N_{i}=S_{i}+I_{i}\) for both females and males so the SIS model can be reduced to
\[ \left \{ \begin{array}{l} \frac{dI_f}{dt} = \beta_{fm} (N_f - I_f(t)) I_m(t)- v_f I_f(t),\\[2ex] \frac{dI_m}{dt} = \beta_{mf} (N_m - I_m(t)) I_f(t)- v_m I_m(t). \end{array} \right. \]
The Gonorrhea model in R
Definition of the model’s papameters
For transmission rates \(\beta_fm=0.000003\) and \(\beta_mf=0.000006\) we have
\[ C= \left( \begin{array}{cc} 0 & 0.000003\\ 0.000006 & 0 \end{array} \right). \]
Note that we assume that for heterosexual intercourse between infected and susceptible individuals, a new infection is twice as likely if the male is the infective. We also assume that the recovery time for females is much longer than the recovery time for males.
parameters <- c(beta1=0.000003,beta2=0.000006,v1=0.007,v2=0.05,N1=10000,N2=15000)
parameters## beta1 beta2 v1 v2 N1 N2
## 3.0e-06 6.0e-06 7.0e-03 5.0e-02 1.0e+04 1.5e+04The state vector
At \(t=0\) there is one infected female and non of the male are infected, i.e., \(I_{f}(0)=1,I_{m}(0)=0\).
state <- c(Y1=1,Y2=0)
state## Y1 Y2
## 1 0times<-seq(0,3000,by=0.1)The transmission model
The R function Gonorrhea is used to implement the transmission model for Gonorrhea in R.
Gonorrhea<-function(t,state,parameters)
{
with(as.list(c(state, parameters)),
{
dY1 <- beta1*(N1-Y1)*Y2-v1*Y1
dY2 <- beta2*(N2-Y2)*Y1-v2*Y2
list(c(dY1,dY2))
})
}Running the model
We use the R function ode() to run the model, the object Gonorrhea.
require(deSolve)
out <- as.data.frame(ode(y=state,times=times,func=Gonorrhea,parms=parameters))
head(out)## time Y1 Y2
## 1 0.0 1.0000000 0.000000000
## 2 0.1 0.9993137 0.008974473
## 3 0.2 0.9986547 0.017898141
## 4 0.3 0.9980229 0.026771503
## 5 0.4 0.9974180 0.035594990
## 6 0.5 0.9968400 0.044369176out1<-out
llast<-length(times)
out1[llast,]## time Y1 Y2
## 30001 3000 7532.051 7121.212Graphical output
Denoting the respective relative removal rates by \(\rho_f = v_f/\beta_f,\) and \(\rho_m = v_m/\beta_m\), the equilibrium values of the above system equal \[ I_f(\infty) = \frac{N_fN_m - \rho_f\rho_m}{\rho_f+N_m} \;\;\; \mbox{and}\;\;\; I_m(\infty) = \frac{N_fN_m - \rho_f\rho_m}{\rho_m+N_f}, \] for females and males, respectively. For our example \(R_{0,f}=4.29\) and \(R_{0,m}=1.80\), The threshold value for transmission should satisfy \(N_{f}N_{m}-\rho_{f}\rho_{m} >0\) or \[ \left ( \frac{N_{f}\beta_{f}}{v_{f}} \times \frac{N_{m}\beta_{m}}{v_{m}} \right ) = R_{0,f} \times R_{0,m} > 1. \] For the above example the basic reproductive number for females and males are \(R_f = 4.29\) and \(R_m = 1.80\), respectively. Note that the product \(R_{0,f} \times R_{0,m}\) should be greater than \(1\).
\[ \begin{array}{l} I_{f}(\infty)=7532.051 \;\;\;\;\mbox{and}\;\;\;\; R_{f}=0.4.29, \\ I_{m}(\infty)=7121.212 \;\;\;\; \mbox{and}\;\;\;\; R_{m}=1.80, \end{array} \] and \(R_{0,f} \times R_{0,m}=7.72\). Figure 1 shows the number of infected individuals in the two populations (male and female) for this scenario.
par(mfrow=c(1,1))
plot(times,out$Y1 ,type="l",main=" ", xlab="time",ylab="-",
ylim=c(0,max(c(out$Y1,out$Y2))))
lines(times,out$Y2,lty=2)
legend(2000,2000,c("Male","Female"),lty=c(1,2))
Figure 1: The Gonorrhea model: number of infected individuals.
Changing the model parameters
Let assume that \(v_{f}\) increases to \(0.0007 \times 5\) while the other parameters in the model remain the same. In this case \[ \begin{array}{l} I_{f}(\infty)=1979.17 \;\;\;\;\mbox{and}\;\;\;\; R_{f}=0.86, \\ I_{m}(\infty)=2878.79 \;\;\;\; \mbox{and}\;\;\;\; R_{m}=1.80, \end{array} \] and \(R_{0,f} \times R_{0,m}=1.54\). In R, the new parameters vector is given by
parameters <- c(beta1=0.000003,beta2=0.000006,v1=0.007*5,v2=0.05,N1=10000,N2=15000)
parameters## beta1 beta2 v1 v2 N1 N2
## 3.0e-06 6.0e-06 3.5e-02 5.0e-02 1.0e+04 1.5e+04Model’s output
We produce the solution for Gonorrhea model with the new parameter setting. Note that, although only the parameter \(v_{f}\) was changed (in the female population) the number of infected individuals was changed also in the male population.
require(deSolve)
out <- as.data.frame(ode(y=state,times=times,func=Gonorrhea,parms=parameters))
#head(out)
out1a<-out
out1a[llast,]## time Y1 Y2
## 30001 3000 1979.167 2878.788Figure 2 shows the number of infected individuals in the two populations for the new scenario.
par(mfrow=c(1,1))
plot(times,out1a$Y1 ,type="l",main=" ", xlab="time",ylab="-",
ylim=c(0,max(c(out1a$Y1,out1a$Y2))))
lines(times,out1a$Y2,lty=2)
legend(1500,1000,c("Male","Female"),lty=c(1,2))
Figure 2: The Gonorrhea model: number of infected individuals.
Changing transmission rates
A new contact matrix
Next, we investigate the influence of \(\beta_{fm}\) on the model. We run the model with a different mixing matrix in which the parameter \(\beta_{mf}\) is fixed and equal to \(\beta_mf=0.000006\) and the parameter \(\beta_fm\) increases and equal to $0.000003 $. The new contact matrix is given by
\[ C= \left( \begin{array}{cc} 0 & 0.000003 \times 5\\ 0.000006 & 0 \end{array} \right). \]
in R, we updated the parameters vector
parameters <- c(beta1=0.000003*5,beta2=0.000006,v1=0.007,v2=0.05,N1=10000,N2=15000)
parameters## beta1 beta2 v1 v2 N1 N2
## 1.5e-05 6.0e-06 7.0e-03 5.0e-02 1.0e+04 1.5e+04while the state vector remains the same.
state <- c(Y1=1,Y2=0)
times<-seq(0,3000,by=0.1)Model’s output
The panels below hows the equilibrium values under the two senarions. Note that,as before, even though that change in the population was related to the female population, the equilibrium value for the male population was changed due to the positive feedback between the female and male populations. This can also be seen in Figure 3 and 4.
require(deSolve)
out2 <- as.data.frame(ode(y=state,times=times,func=Gonorrhea,parms=parameters))
out1[llast,]## time Y1 Y2
## 30001 3000 7532.051 7121.212out2[llast,]## time Y1 Y2
## 30001 3000 9446.839 7969.697par(mfrow=c(1,1))
plot(times,out1$Y1 ,type="l",main="Female", xlab="time",ylab="-",
ylim=c(0,max(c(out1$Y1,out2$Y1))))
lines(times,out2$Y1,lty=2)
legend(1000,2000,c("beta_fm=0.000003","beta_fm=0.000003*5"),lty=c(1,2))
Figure 3: The Gonorrhea model: female population.
par(mfrow=c(1,1))
plot(times,out1$Y2 ,type="l",main="Male", xlab="time",ylab="-",
ylim=c(0,max(c(out1$Y2,out2$Y2))))
lines(times,out2$Y2,lty=2)
legend(1000,2000,c("beta_fm=0.000003","beta_fm=0.000003*5"),lty=c(1,2))
Figure 4: The Gonorrhea model: male population.
3. Example 2: SIR Model With Two Interacting Sub-populations
The mixing matrix
The Gonorrhea model discussed in the previous section assumes that there is no mixing within each sub-population (\(\beta_{mm}=\beta_{ff}=0\)). In this section we consider a SIR model in which individuals are assumed to make contact within and across sub-populations. For a case with two sub-populations the WAIFW or mixing matrix is given by
\[ C= \left( \begin{array}{cc} \beta_{11} & \beta_{12}\\ \beta_{21} & \beta_{22} \end{array} \right). \]
For each sub-population, the population size is \(N_{i}=S_{i}+I_{i}+R_{i}\), \(i=1,2\).
The number of new cases and the force of infection
The main difference between the current model to the SIR model for a single population is that in the current model the number of new cases depends on the mixing pattern within and across sub-populations, that is \[ \left( \Sigma_{j=1}^{2}\beta_{ij}I_{j} \right )S_{i}=\lambda_{i} S_{i}. \]
The SIR model
The differential equation system for the i\(th\) subpopulation is given by \[ \left \{ \begin{array}{l} \frac{dS_{i}(t)}{dt} =- \left( \Sigma_{j=1}^{2}\beta_{ij}I_{j} \right )S_{i} +N_{i}\mu_{i}-\mu_{i}S_{i} ,\\ \frac{dI_{i}(t)}{dt} = \left( \Sigma_{j=1}^{2}\beta_{ij}I_{j} \right )S_{i} - (\mu_{i}+\nu_{i}) I_{i} ,\\ \frac{dR_{i}(t)}{dt} = \nu_{i}I_{i}-\mu_{i}R_{i}. \end{array} \right. \]
Implementation in R
Let us consider the special in which the mixing matrix in the population given by
\[ C= \left( \begin{array}{cc} 0.05 & 0.05\\ 0.075 & 0.05 \end{array} \right). \]
Since there are two populations we need to specify 6 differential equations in the function SIRtwo. The R objects Y1, Y2 and Y3 correspond to S, I and R in the first population and Y4, Y5 and Y6 correspond to S, I and R in the second population.
The SIRtwo function
SIRtwo<-function(t,state,parameters)
{
with(as.list(c(state, parameters)),
{
dY1 <- -(beta11*Y2+beta12*Y5)*Y1+mu-mu*Y1
dY2 <- (beta11*Y2+beta12*Y5)*Y1-v1*Y2-mu*Y2
dY3 <- v1*Y2 - mu*Y3
dY4 <- -(beta21*Y2+beta22*Y5)*Y4+mu-mu*Y4
dY5 <- (beta21*Y2+beta22*Y5)*Y4-v2*Y5-mu*Y5
dY6 <- v2*Y5-mu*Y6
list(c(dY1,dY2,dY3,dY4,dY5,dY6))
})
}Model parameters and state variables
We assume that at \(t=0\), \(80\%\) of the individuals in the two populations are susceptible \(S_{1}(0)=S_{2}(0)=0.8\) and \(20\%\) are infectious, \(I_{1}(0)=I_{2}(0)=0.2\),
state <- c(Y1=0.8,Y2=0.2,Y3=0,Y4=0.8,Y5=0.2,Y6=0)
state## Y1 Y2 Y3 Y4 Y5 Y6
## 0.8 0.2 0.0 0.8 0.2 0.0times<-seq(0,4000,by=0.01)For \(\mu=0.001\) and \(\nu_{1}=\nu_{2}=\frac{1}{30}\), the parameters vector in R is given by
#parameters <- c(beta11=0.05,beta12=0.075,beta21=0.075,beta22=0.05,v1=1/30,v2=1/30,mu=0.001) #C: symetric
parameters <- c(beta11=0.05,beta12=0.05,beta21=0.075,beta22=0.05,v1=1/30,v2=1/30,mu=0.001) #basic setting
#parameters <- c(beta11=0.05,beta12=0.0,beta21=0.075,beta22=0.05,v1=1/30,v2=1/30,mu=0.001) #beta12=0.00
parameters ## beta11 beta12 beta21 beta22 v1 v2 mu
## 0.05000000 0.05000000 0.07500000 0.05000000 0.03333333 0.03333333 0.00100000Running the model
require(deSolve)
out <- as.data.frame(ode(y=state,times=times,func=SIRtwo,parms=parameters))
head(out)## time Y1 Y2 Y3 Y4 Y5 Y6
## 1 0.00 0.8000000 0.2000000 0.000000e+00 0.8000000 0.2000000 0.000000e+00
## 2 0.01 0.7998420 0.2000913 6.668156e-05 0.7998020 0.2001313 6.668822e-05
## 3 0.02 0.7996839 0.2001827 1.333929e-04 0.7996039 0.2002627 1.334196e-04
## 4 0.03 0.7995257 0.2002741 2.001340e-04 0.7994057 0.2003941 2.001940e-04
## 5 0.04 0.7993676 0.2003655 2.669050e-04 0.7992076 0.2005254 2.670116e-04
## 6 0.05 0.7992093 0.2004570 3.337057e-04 0.7990093 0.2006568 3.338723e-04out1<-outGraphical output (\(\beta_{12}=0\))
Figure 5 shows the number of susceptible in the two populations.
par(mfrow=c(2,1))
plot(times,out1$Y1,type="l")
title("Population 1")
plot(times,out1$Y4,type="l")
title("population 2")
Figure 5: SIR model with two interacting populations.
Changing the transmiison rates (I)
The new contact matrix with \(\beta_{12}=0.075\)
In this example, \(\beta_{12}=0.075\). Hence, the contact metrix is symetric and given by
\[ C= \left( \begin{array}{cc} 0.05 & 0.075\\ 0.075 & 0.05 \end{array} \right). \]
In R, the new parameters vector is given by
parameters <- c(beta11=0.05,beta12=0.075,beta21=0.075,beta22=0.05,v1=1/30,v2=1/30,mu=0.001) #basic setting
parameters## beta11 beta12 beta21 beta22 v1 v2 mu
## 0.05000000 0.07500000 0.07500000 0.05000000 0.03333333 0.03333333 0.00100000Running the model
require(deSolve)
out <- as.data.frame(ode(y=state,times=times,func=SIRtwo,parms=parameters))
head(out)## time Y1 Y2 Y3 Y4 Y5 Y6
## 1 0.00 0.8000000 0.2000000 0.000000e+00 0.8000000 0.2000000 0.000000e+00
## 2 0.01 0.7998020 0.2001314 6.668823e-05 0.7998020 0.2001314 6.668823e-05
## 3 0.02 0.7996038 0.2002627 1.334196e-04 0.7996038 0.2002627 1.334196e-04
## 4 0.03 0.7994056 0.2003942 2.001941e-04 0.7994056 0.2003942 2.001941e-04
## 5 0.04 0.7992074 0.2005256 2.670117e-04 0.7992074 0.2005256 2.670117e-04
## 6 0.05 0.7990090 0.2006571 3.338725e-04 0.7990090 0.2006571 3.338725e-04out2<-outModel’s output
Figure 6 and 7 show the number of susceptible under the two scenarios in the two populations. Note that the inter-epidemic period under scenario 2 (\(\beta=0.0075\) is shooter due to the higher transmission rate from population 2 to 1.
par(mfrow=c(1,1))
plot(times,out1$Y1,type="l",ylim=c(0,1))
lines(times,out2$Y1,lty=2,col=2)
legend(0,0.95,c("beta_12=0.075","beta_12=0.05"),lty=c(2,1),col=c(2,1))
title("susceptible in population 1")
Figure 6: SIR model with two interacting populations: susceptible in population 1.
par(mfrow=c(1,1))
plot(times,out1$Y4,type="l",ylim=c(0,1))
lines(times,out2$Y4,lty=2,col=2)
legend(0,0.95,c("beta_12=0.075","beta_12=0.05"),lty=c(2,1),col=c(2,1))
title("susceptible in population 2")
Figure 7: SIR model with two interacting populations: susceptible in population 2.
Changing the transmiison rates (II)
The new contact matrix with \(\beta_{12}=0.0\)
In thi example, \(\beta_{12}=0\). Hence, the contact metrix is given by
\[ C= \left( \begin{array}{cc} 0.05 & 0.0\\ 0.075 & 0.05 \end{array} \right). \]
and the new parameters vector
parameters <- c(beta11=0.05,beta12=0.0,beta21=0.075,beta22=0.05,v1=1/30,v2=1/30,mu=0.001) #basic setting
parameters## beta11 beta12 beta21 beta22 v1 v2 mu
## 0.05000000 0.00000000 0.07500000 0.05000000 0.03333333 0.03333333 0.00100000Running the model
require(deSolve)
out <- as.data.frame(ode(y=state,times=times,func=SIRtwo,parms=parameters))
head(out)## time Y1 Y2 Y3 Y4 Y5 Y6
## 1 0.00 0.8000000 0.2000000 0.000000e+00 0.8000000 0.2000000 0.000000e+00
## 2 0.01 0.7999220 0.2000113 6.666822e-05 0.7998020 0.2001313 6.668822e-05
## 3 0.02 0.7998440 0.2000227 1.333396e-04 0.7996040 0.2002626 1.334195e-04
## 4 0.03 0.7997660 0.2000340 2.000140e-04 0.7994060 0.2003938 2.001939e-04
## 5 0.04 0.7996880 0.2000453 2.666915e-04 0.7992079 0.2005251 2.670114e-04
## 6 0.05 0.7996101 0.2000566 3.333722e-04 0.7990099 0.2006562 3.338720e-04out3<-outModel’s output
Figure 8 and 9 show that, as expected, since the transmission rate from population 1 to 2 is equal to 0, the inter-epidemic period is longer under the new scenario.
par(mfrow=c(1,1))
plot(times,out3$Y1,type="l",ylim=c(0,1))
lines(times,out2$Y1,lty=2,col=2)
legend(0,0.95,c("beta_12=0.075","beta_12=0.0"),lty=c(2,1),col=c(2,1))
title("susceptible in population 1")
Figure 8: SIR model with two interacting populations: susceptible in population 1.
par(mfrow=c(1,1))
plot(times,out3$Y4,type="l",ylim=c(0,1))
lines(times,out2$Y4,lty=2,col=2)
legend(0,0.95,c("beta_12=0.075","beta_12=0.0"),lty=c(2,1),col=c(2,1))
title("susceptible in population 2")
Figure 9: SIR model with two interacting populations: susceptible in population 2.
Changing the transmiison rates (III)
The new contact matrix with \(\beta_{12}=0.0\)
In this example, there is not transmissionm for a contact of two individuals from population 1, \(\beta_{11}=0\). Hence, the contact metrix is given by
\[ C= \left( \begin{array}{cc} 0.0 & 0.007\\ 0.075 & 0.05 \end{array} \right). \]
and the parameters vector by
parameters <- c(beta11=0.0,beta12=0.075,beta21=0.075,beta22=0.05,v1=1/30,v2=1/30,mu=0.001) #basic setting
parameters## beta11 beta12 beta21 beta22 v1 v2 mu
## 0.00000000 0.07500000 0.07500000 0.05000000 0.03333333 0.03333333 0.00100000Running the model
require(deSolve)
out <- as.data.frame(ode(y=state,times=times,func=SIRtwo,parms=parameters))
head(out)## time Y1 Y2 Y3 Y4 Y5 Y6
## 1 0.00 0.8000000 0.2000000 0.000000e+00 0.8000000 0.2000000 0.000000e+00
## 2 0.01 0.7998820 0.2000514 6.667489e-05 0.7998020 0.2001313 6.668822e-05
## 3 0.02 0.7997639 0.2001028 1.333662e-04 0.7996039 0.2002626 1.334195e-04
## 4 0.03 0.7996457 0.2001542 2.000741e-04 0.7994059 0.2003939 2.001940e-04
## 5 0.04 0.7995275 0.2002057 2.667984e-04 0.7992077 0.2005252 2.670115e-04
## 6 0.05 0.7994093 0.2002572 3.335392e-04 0.7990096 0.2006565 3.338721e-04out4<-outModel’s output
Figure 10 and 11 show that the inter-epidemic period is longer in both populations even though, the change is only in the transmission is population 1.
par(mfrow=c(1,1))
plot(times,out4$Y1,type="l",ylim=c(0,1))
lines(times,out2$Y1,lty=2,col=2)
legend(0,0.95,c("beta_11=0.05","beta_11=0.0"),lty=c(2,1),col=c(2,1))
title("susceptible in population 1")
Figure 10: SIR model with two interacting populations: susceptible in population 1.
par(mfrow=c(1,1))
plot(times,out4$Y4,type="l",ylim=c(0,1))
lines(times,out2$Y4,lty=2,col=2)
legend(0,0.95,c("beta_11=0.05","beta_11=0.0"),lty=c(2,1),col=c(2,1))
title("susceptible in population 2")
Figure 11: SIR model with two interacting populations: susceptible in population 2.