18 February 2022
Julien Arino
Department of Mathematics & Data Science Nexus University of Manitoba*
Canadian Centre for Disease Modelling Canadian COVID-19 Mathematical Modelling Task Force NSERC-PHAC EID Modelling Consortium (CANMOD, MfPH, OMNI/RÉUNIS)
It first began, it is said, in the parts of Ethiopia above Egypt, and thence descended into Egypt and Libya and into most of the King's country [Persia]. Suddenly falling upon Athens, it first attacked the population in Piraeus—which was the occasion of their saying that the Peloponnesians had poisoned the reservoirs, there being as yet no wells there—and afterwards appeared in the upper city, when the deaths became much more frequent.
History of the Peloponnesian War
Mobility is complicated:
And yet mobility drives spatio-temporal spread:
Pathogens travel along trade routes
In ancient times, trade routes were relatively easy to comprehend
With acceleration and globalization of mobility, things change
All migrants/travellers carry with them their "health history"
Pathogens ignore borders and politics
Index case for international spread arrives HKG 21 February 2003
Last country with local transmission (Taiwan) removed from list 5 July 2003
8273 cases in 28 countries
(Of these cases, 1706 were HCW)
775 deaths (CFR 9.4%)
JA. Describing, modelling and forecasting the spatial and temporal spread of COVID-19 - A short review. Fields Institute Communications 85:25-51 (2022)
By 29 January, virus was found in all provinces of mainland China
⟹ \implies⟹ amplification must have been occurring for a while longer
Disease spread process in a jurisdiction-based world
NcpN_{cp}Ncp population of compartment c∈Cc\in\mathcal{C}c∈C in patch p∈Pp\in\mathcal{P}p∈P
Assume no birth or death. Balance inflow and outflow
Ncp′=(∑q∈P∖{p}mcpqNcq)−(∑q∈P∖{p}mcqp)Ncp=∑q∈PmcpqNcqwith mcpp=−∑q∈P∖{p}mcqp\begin{aligned} N_{cp}' &= \left(\sum_{q\in\mathcal{P}\setminus\{p\}} m_{cpq}N_{cq}\right)-\left(\sum_{q\in\mathcal{P}\setminus\{p\}} m_{cqp}\right)N_{cp} \\ &\\ &= \sum_{q\in\mathcal{P}} m_{cpq}N_{cq} \qquad\textrm{with } m_{cpp}=-\sum_{q\in\mathcal{P}\setminus\{p\}} m_{cqp} \end{aligned} Ncp′=⎝⎛q∈P∖{p}∑mcpqNcq⎠⎞−⎝⎛q∈P∖{p}∑mcqp⎠⎞Ncp=q∈P∑mcpqNcqwith mcpp=−q∈P∖{p}∑mcqp
B(N)B(N)B(N) is the birth rate (typically bbb or bNbNbN)
LLL = latently infected (≃E\simeq E≃E exposed, although the latter term is ambiguous)
Sp′=Bp(Np)+νpRp−Φp−dpSp+∑q∈PmSpqSqLp′=Φp−(εp+dp)Lp+∑q∈PmLpqLqIp′=εpLp−(γp+dp)Ip+∑q∈PmIpqIqRp′=γpIp−(νp+dp)Rp+∑q∈PmRpqRq\begin{aligned} S_{p}' &=\mathcal{B}_p\left(N_p\right)+\nu_pR_p-\Phi_p-d_pS_p \red{+\textstyle{\sum_{q\in\mathcal{P}}} m_{Spq}S_{q}} \\ L_{p}' &=\Phi_p-\left( \varepsilon_{p}+d_{p}\right)L_{p} \red{+\textstyle{\sum_{q\in\mathcal{P}}} m_{Lpq}L_{q}} \\ I_{p}' &=\varepsilon_pL_p-(\gamma_p+d_p)I_p \red{+\textstyle{\sum_{q\in\mathcal{P}}} m_{Ipq}I_{q}} \\ R_{p}' &=\gamma _{p}I_{p}-\left(\nu_{p}+d_{p}\right)R_{p} \red{+\textstyle{\sum_{q\in\mathcal{P}}} m_{Rpq}R_{q}} \end{aligned} Sp′Lp′Ip′Rp′=Bp(Np)+νpRp−Φp−dpSp+∑q∈PmSpqSq=Φp−(εp+dp)Lp+∑q∈PmLpqLq=εpLp−(γp+dp)Ip+∑q∈PmIpqIq=γpIp−(νp+dp)Rp+∑q∈PmRpqRq
Φp=βpSpIpNpqp,qp∈{0,1}\Phi_p=\beta_p\frac{S_pI_p}{N_p^{q_p}},\qquad q_p\in\{0,1\} Φp=βpNpqpSpIp,qp∈{0,1}
S\mathcal{S}S a set of species
Ssp′=Bsp(Nsp)+νspRsp−Φsp−dspSsp+∑q∈PmSspqSsqLsp′=Φsp−(εsp+dsp)Lsp+∑q∈PmLspqLsqIsp′=εspLsp−(γsp+dsp)Isp+∑q∈PmIspqIsqRsp=γspIsp−(νsp+dsp)Rsp+∑q∈PmRspqRsq\begin{aligned} S_{sp}' &= \mathcal{B}_{sp}(N_{sp})+\nu_{sp}R_{sp}-\Phi_{sp}-d_{sp}S_{sp} \red{+\textstyle{\sum_{q\in\mathcal{P}}} m_{Sspq}S_{sq}} \\ L_{sp}' &= \Phi_{sp}-(\varepsilon_{sp}+d_{sp})L_{sp} \red{+\textstyle{\sum_{q\in\mathcal{P}}}m_{Lspq}L_{sq}} \\ I_{sp}' &= \varepsilon_{sp}L_{sp}-(\gamma_{sp}+d_{sp})I_{sp} \red{+\textstyle{\sum_{q\in\mathcal{P}}} m_{Ispq}I_{sq}} \\ R_{sp} &= \gamma _{sp}I_{sp}-(\nu_{sp}+d_{sp})R_{sp} \red{+\textstyle{\sum_{q\in\mathcal{P}}} m_{Rspq}R_{sq}} \end{aligned} Ssp′Lsp′Isp′Rsp=Bsp(Nsp)+νspRsp−Φsp−dspSsp+∑q∈PmSspqSsq=Φsp−(εsp+dsp)Lsp+∑q∈PmLspqLsq=εspLsp−(γsp+dsp)Isp+∑q∈PmIspqIsq=γspIsp−(νsp+dsp)Rsp+∑q∈PmRspqRsq
Φsp=∑k∈SβskpSspIkpNpqp,qp∈{0,1}\Phi_{sp}=\sum_{k\in\mathcal{S}}\beta_{skp}\frac{S_{sp}I_{kp}}{N_p^{q_p}},\qquad q_p\in\{0,1\} Φsp=k∈S∑βskpNpqpSspIkp,qp∈{0,1}
Spq′=Bpq(Npr)+νpqRpq−Φpq−dpqSpq+∑k∈PmSpqkSpkLpq′=Φpq−(εpq+dpq)Lpq+∑k∈PmLpqkLpkIpq′=εpqLpq−(γpq+dpq)Ipq+∑k∈PmIpqkIpkRpq′=γpqIpq−(νpq+dpq)Rpq+∑k∈PmRpqkRpk\begin{aligned} S_{pq}' =& \mathcal{B}_{pq}\left(N_p^r\right)+\nu_{pq} R_{pq}-\Phi_{pq}-d_{pq}S_{pq} \red{+\textstyle{\sum_{k\in\mathcal{P}}} m_{Spqk}S_{pk}} \\ L_{pq}' =& \Phi_{pq} -(\varepsilon_{pq}+d_{pq})L_{pq} \red{+\textstyle{\sum_{k\in\mathcal{P}}} m_{Lpqk}L_{pk}} \\ I_{pq}' =& \varepsilon_{pq} L_{pq} -(\gamma_{pq}+d_{pq})I_{pq} \red{+\textstyle{\sum_{k\in\mathcal{P}}} m_{Ipqk}I_{pk}} \\ R_{pq}' =& \gamma_{pq} I_{pq} -(\nu_{pq}+d_{pq})R_{pq} \red{+\textstyle{\sum_{k\in\mathcal{P}}} m_{Rpqk}R_{pk}} \end{aligned} Spq′=Lpq′=Ipq′=Rpq′=Bpq(Npr)+νpqRpq−Φpq−dpqSpq+∑k∈PmSpqkSpkΦpq−(εpq+dpq)Lpq+∑k∈PmLpqkLpkεpqLpq−(γpq+dpq)Ipq+∑k∈PmIpqkIpkγpqIpq−(νpq+dpq)Rpq+∑k∈PmRpqkRpk
Φpq=∑k∈PβpqkSpqIkqNpqq,qq={0,1}\Phi_{pq}=\sum_{k\in\mathcal{P}}\beta_{pqk}\frac{S_{pq}I_{kq}}{N_p^{q_q}},\qquad q_q=\{0,1\} Φpq=k∈P∑βpqkNpqqSpqIkq,qq={0,1}
U⊊C\mathcal{U}\subsetneq\mathcal{C}U⊊C uninfected and I⊊C\mathcal{I}\subsetneq\mathcal{C}I⊊C infected compartments, U∪I=C\mathcal{U}\cup\mathcal{I}=\mathcal{C}U∪I=C and U∩I=∅\mathcal{U}\cap\mathcal{I}=\emptysetU∩I=∅
For k∈Uk\in\mathcal{U}k∈U, ℓ∈I\ell\in\mathcal{I}ℓ∈I and p∈Pp\in\mathcal{P}p∈P,
skp′=fkp(Sp,Ip)+∑q∈Pmkpqskqiℓp′=gℓp(Sp,Ip)+∑q∈Pmℓpqiℓq\begin{aligned} s_{kp}' &= f_{kp}(S_p,I_p)+\sum_{q\in\mathcal{P}} m_{kpq}s_{kq} \\ i_{\ell p}' &= g_{\ell p}(S_p,I_p)+\sum_{q\in\mathcal{P}} m_{\ell pq}i_{\ell q} \end{aligned} skp′iℓp′=fkp(Sp,Ip)+q∈P∑mkpqskq=gℓp(Sp,Ip)+q∈P∑mℓpqiℓq
where Sp=(s1p,…,s∣U∣p)S_p=(s_{1p},\ldots,s_{|\mathcal{U}|p})Sp=(s1p,…,s∣U∣p) and Ip=(i1p,…,i∣I∣p)I_p=(i_{1p},\ldots,i_{|\mathcal{I}|p})Ip=(i1p,…,i∣I∣p)
For simplicity, consider ∣P∣|\mathcal{P}|∣P∣-SLIRS with Bp(Np)=Bp\mathcal{B}_p(N_p)=\mathcal{B}_pBp(Np)=Bp
Sp′=Bp−Φp−dpSp+νpRp+∑q∈PmSpqSqLp′=Φp−(εp+dp)Lp+∑q∈PmLpqLqIp′=εpLp−(γp+dp)Ip+∑q∈PmIpqIqRp′=γpIp−(νp+dp)Rp+∑q∈PmRpqRq\begin{aligned} S_{p}' &=\mathcal{B}_p-\Phi_p-d_pS_p+\nu_pR_p +\textstyle{\sum_{q\in\mathcal{P}}} m_{Spq}S_{q} \\ L_{p}' &=\Phi_p-\left( \varepsilon_{p}+d_{p}\right)L_{p} +\textstyle{\sum_{q\in\mathcal{P}}} m_{Lpq}L_{q} \\ I_{p}' &=\varepsilon_pL_p-(\gamma_p+d_p)I_p +\textstyle{\sum_{q\in\mathcal{P}}} m_{Ipq}I_{q} \\ R_{p}' &=\gamma _{p}I_{p}-\left(\nu_{p}+d_{p}\right)R_{p} +\textstyle{\sum_{q\in\mathcal{P}}} m_{Rpq}R_{q} \end{aligned} Sp′Lp′Ip′Rp′=Bp−Φp−dpSp+νpRp+∑q∈PmSpqSq=Φp−(εp+dp)Lp+∑q∈PmLpqLq=εpLp−(γp+dp)Ip+∑q∈PmIpqIq=γpIp−(νp+dp)Rp+∑q∈PmRpqRq
System of 4∣P∣4|\mathcal{P}|4∣P∣ equations
System of 4∣P∣4|\mathcal{P}|4∣P∣ equations !!!
However, a lot of structure:
⟹ \implies⟹ Good case of large-scale system (matrix analysis is your friend)
M∈Mn(R)=Rn×nM\in\mathcal{M}_n(\mathbb{R})=\mathbb{R}^{n\times n}M∈Mn(R)=Rn×n a square matrix with entries denoted mijm_{ij}mij
M≥0M\geq\mathbf{0}M≥0 if mij≥0m_{ij}\geq 0mij≥0 for all i,ji,ji,j (could be the zero matrix); M>0M>\mathbf{0}M>0 if M≥0M\geq\mathbf{0}M≥0 and ∃i,j\exists i,j∃i,j with mij>0m_{ij}>0mij>0; M≫0M\gg\mathbf{0}M≫0 if mij>0m_{ij}>0mij>0 ∀i,j=1,…,n\forall i,j=1,\ldots,n∀i,j=1,…,n. Same notation for vectors
σ(M)={λ∈C;Mλ=λv,v≠0}\sigma(M)=\{\lambda\in\mathbf{C}; M\lambda=\lambda\mathbf{v}, \mathbf{v}\neq\mathbf{0}\}σ(M)={λ∈C;Mλ=λv,v=0} spectrum of MMM
ρ(M)=maxλ∈σ(M){∣λ∣}\rho(M)=\max_{\lambda\in\sigma(M)}\{|\lambda|\}ρ(M)=maxλ∈σ(M){∣λ∣} spectral radius
s(M)=maxλ∈σ(M){ℜ(λ)}s(M)=\max_{\lambda\in\sigma(M)}\{\Re(\lambda)\}s(M)=maxλ∈σ(M){ℜ(λ)} spectral abscissa (or stability modulus)
MMM is an M-matrix if it is a Z-matrix (mij≤0m_{ij}\leq 0mij≤0 for i≠ji\neq ji=j) and M=sI−AM = s\mathbb{I}-AM=sI−A, with A≥0A\geq 0A≥0 and s≥ρ(A)s\geq \rho(A)s≥ρ(A)
Consider behaviour of Np=Sp+Lp+Ip+RpN_p=S_p+L_p+I_p+R_pNp=Sp+Lp+Ip+Rp. We have
Np′=Bp−Φp−dpSp+νpRp+∑q∈PmSpqSq+Φp−(εp+dp)Lp+∑q∈PmLpqLq+εpLp−(γp+dp)Ip+∑q∈PmIpqIq+γpIp−(νp+dp)Rp+∑q∈PmRpqRq\begin{aligned} N_p' &=\mathcal{B}_p\cancel{-\Phi_p}-d_pS_p\cancel{+\nu_pR_p} +\textstyle{\sum_{q\in\mathcal{P}}} m_{Spq}S_{q} \\ &\quad \cancel{+\Phi_p}-\left(\cancel{\varepsilon_{p}} +d_{p}\right)L_{p} +\textstyle{\sum_{q\in\mathcal{P}}} m_{Lpq}L_{q} \\ &\quad \cancel{+\varepsilon_pL_p}-(\cancel{\gamma_p}+d_p)I_p +\textstyle{\sum_{q\in\mathcal{P}}} m_{Ipq}I_{q} \\ &\quad \cancel{+\gamma _{p}I_{p}} -\left(\cancel{\nu_{p}}+d_{p}\right)R_{p} +\textstyle{\sum_{q\in\mathcal{P}}} m_{Rpq}R_{q} \end{aligned} Np′=Bp−Φp−dpSp+νpRp+∑q∈PmSpqSq+Φp−(εp+dp)Lp+∑q∈PmLpqLq+εpLp−(γp+dp)Ip+∑q∈PmIpqIq+γpIp−(νp+dp)Rp+∑q∈PmRpqRq
So
Np′=Bp−dpNp+∑X∈{S,L,I,R}∑q∈PmXpqXqN_p'=\mathcal{B}_p-d_pN_p +\sum_{X\in\{S,L,I,R\}}\sum_{q\in\mathcal{P}} m_{Xpq}X_{q} Np′=Bp−dpNp+X∈{S,L,I,R}∑q∈P∑mXpqXq
We have
Write this in vector form
N′=b−dN+∑X∈{S,L,I,R}MXX\mathbf{N}'=\mathbf{b}-\mathbf{d}\mathbf{N}+\sum_{X\in\{S,L,I,R\}}\mathcal{M}^X\mathbf{X} N′=b−dN+X∈{S,L,I,R}∑MXX
where b=(B1,…,B∣P∣)T,N=(N1,…,N∣P∣)T,X=(X1,…,X∣P∣)T∈R∣P∣,\mathbf{b}=(\mathcal{B}_1,\ldots,\mathcal{B}_{|\mathcal{P}|})^T,\mathbf{N}=(N_1,\ldots,N_{|\mathcal{P}|})^T,\mathbf{X}=(X_1,\ldots,X_{|\mathcal{P}|})^T\in\mathbb{R}^{|\mathcal{P}|},b=(B1,…,B∣P∣)T,N=(N1,…,N∣P∣)T,X=(X1,…,X∣P∣)T∈R∣P∣, d,MX\mathbf{d},\mathcal{M}^Xd,MX ∣P∣×∣P∣|\mathcal{P}|\times|\mathcal{P}|∣P∣×∣P∣-matrices with
d=diag(d1,…,d∣P∣)\mathbf{d}=\mathsf{diag}\left(d_1,\ldots,d_{|\mathcal{P}|}\right) d=diag(d1,…,d∣P∣)
Mc=(−∑q∈Pmcq1mc12mc1∣P∣mc21−∑q∈Pmcq2mc2∣P∣mc∣P∣1mc∣P∣2−∑q∈Pmcq∣P∣)\mathcal{M}^c= \begin{pmatrix} -\sum_{q\in\mathcal{P}} m_{cq1} & m_{c12} & & m_{c1|\mathcal{P}|} \\ m_{c21} & -\sum_{q\in\mathcal{P}} m_{cq2} & & m_{c2|\mathcal{P}|} \\ & & & \\ m_{c|\mathcal{P}|1} & m_{c|\mathcal{P}|2} & & -\sum_{q\in\mathcal{P}} m_{cq|\mathcal{P}|} \end{pmatrix} Mc=⎝⎛−∑q∈Pmcq1mc21mc∣P∣1mc12−∑q∈Pmcq2mc∣P∣2mc1∣P∣mc2∣P∣−∑q∈Pmcq∣P∣⎠⎞
Consider a compartment c∈Cc\in\mathcal{C}c∈C. Then the following hold true:
Recall that
Suppose movement rates equal for all compartments, i.e.,
MS=ML=MI=MR=:M\mathcal{M}^S=\mathcal{M}^L=\mathcal{M}^I=\mathcal{M}^R=:\mathcal{M} MS=ML=MI=MR=:M
Then
N′=b−dN+M∑X∈{S,L,I,R}X=b−dN+MN\begin{aligned} \mathbf{N}' &= \mathbf{b}-\mathbf{d}\mathbf{N}+\mathcal{M}\sum_{X\in\{S,L,I,R\}}\mathbf{X}\\ &= \mathbf{b}-\mathbf{d}\mathbf{N}+\mathcal{M}\mathbf{N} \end{aligned} N′=b−dN+MX∈{S,L,I,R}∑X=b−dN+MN
N′=b−dN+MN\mathbf{N}'=\mathbf{b}-\mathbf{d}\mathbf{N}+\mathcal{M}\mathbf{N} N′=b−dN+MN
Equilibria
N′=0⇔b−dN+MN=0⇔(d−M)N=b⇔N⋆=(d−M)−1b\begin{aligned} \mathbf{N}'=\mathbf{0} &\Leftrightarrow \mathbf{b}-\mathbf{d}\mathbf{N}+\mathcal{M}\mathbf{N}=\mathbf{0} \\ &\Leftrightarrow (\mathbf{d}-\mathcal{M})\mathbf{N}=\mathbf{b} \\ &\Leftrightarrow \mathbf{N}^\star=(\mathbf{d}-\mathcal{M})^{-1}\mathbf{b} \end{aligned} N′=0⇔b−dN+MN=0⇔(d−M)N=b⇔N⋆=(d−M)−1b
given, of course, that d−M\mathbf{d}-\mathcal{M}d−M (or, equivalently, M−d\mathcal{M}-\mathbf{d}M−d) is invertible.. Is it?
M\mathcal{M}M a movement matrix and DDD a diagonal matrix. The following hold true:
Using a spectrum shift,
s(M−d)=−minp∈Pdps(\mathcal{M}-\mathbf{d})=-\min_{p\in\mathcal{P}}d_p s(M−d)=−p∈Pmindp
This gives a constraint: for total population to behave well (in general, we want this), we must assume all death rates are positive
Assume they are (in other words, assume d\mathbf{d}d nonsingular). Then M−d\mathcal{M}-\mathbf{d}M−d is nonsingular and N⋆=(d−M)−1b\mathbf{N}^\star=(\mathbf{d}-\mathcal{M})^{-1}\mathbf{b}N⋆=(d−M)−1b unique
N⋆=(d−M)−1b\mathbf{N}^\star=(\mathbf{d}-\mathcal{M})^{-1}\mathbf{b}N⋆=(d−M)−1b attracts solutions of
N′=b−dN+MN=:f(N)\mathbf{N}'=\mathbf{b}-\mathbf{d}\mathbf{N}+\mathcal{M}\mathbf{N}=:f(\mathbf{N}) N′=b−dN+MN=:f(N)
Indeed, we have
Df=M−dDf=\mathcal{M}-\mathbf{d} Df=M−d
Since we now assume that d\mathbf{d}d is nonsingular, we have (spectral shift & properties of M\mathcal{M}M) s(M−d)=−minp∈Pdp<0s(\mathcal{M}-\mathbf{d})=-\min_{p\in\mathcal{P}}d_p<0s(M−d)=−minp∈Pdp<0
M\mathcal{M}M irreducible →\rightarrow→ N⋆≫0\mathbf{N}^\star\gg 0N⋆≫0 (provided b>0\mathbf{b}>\mathbf{0}b>0, of course)
Assume M\mathcal{M}M reducible. Let aaa be the number of minimal absorbing sets in the corresponding connection graph G(M)\mathcal{G}(\mathcal{M})G(M). Then
From Foster and Jacquez, Multiple zeros for eigenvalues and the multiplicity of traps of a linear compartmental system, Mathematical Biosciences (1975)
Suppose movement rates similar for all compartments, i.e., the zero/nonzero patterns in all matrices are the same but not the entries
Let
M‾=[minX∈{S,L,I,R}mXpq]pq,p≠qM‾=[maxX∈{S,L,I,R}mXpq]pq,p=q\underline{\mathcal{M}}=\left[\min_{X\in\{S,L,I,R\}}m_{Xpq}\right]_{pq,p\neq q}\qquad \underline{\mathcal{M}}=\left[\max_{X\in\{S,L,I,R\}}m_{Xpq}\right]_{pq,p=q} M=[X∈{S,L,I,R}minmXpq]pq,p=qM=[X∈{S,L,I,R}maxmXpq]pq,p=q
and
M‾=[maxX∈{S,L,I,R}mXpq]pq,p≠qM‾=[minX∈{S,L,I,R}mXpq]pq,p=q\overline{\mathcal{M}}=\left[\max_{X\in\{S,L,I,R\}}m_{Xpq}\right]_{pq,p\neq q}\qquad \overline{\mathcal{M}}=\left[\min_{X\in\{S,L,I,R\}}m_{Xpq}\right]_{pq,p=q} M=[X∈{S,L,I,R}maxmXpq]pq,p=qM=[X∈{S,L,I,R}minmXpq]pq,p=q
Then we have
b−dN+M‾N≤N′≤b−dN+M‾N\mathbf{b}-\mathbf{d}\mathbf{N}+\underline{\mathcal{M}}\mathbf{N}\leq\mathbf{N}'\leq\mathbf{b}-\mathbf{d}\mathbf{N}+\overline{\mathcal{M}}\mathbf{N} b−dN+MN≤N′≤b−dN+MN
Me, roughly every 6 months: Oooh, coooool, a linear differential inclusion!
Me, roughly 10 minutes after that previous statement: Quel con!
Indeed M‾\underline{\mathcal{M}}M and M‾\overline{\mathcal{M}}M are are not movement matrices (in particular, their column sums are not all zero)
So no luck there..
However, non lasciate ogne speranza, we can still do stuff!
Assume system at equilibrium and Lp=Ip=0L_p=I_p=0Lp=Ip=0 for p∈Pp\in\mathcal{P}p∈P. Then Φp=0\Phi_p=0Φp=0 and
0=Bp−dpSp+νpRp+∑q∈PmSpqSq0=−(νp+dp)Rp+∑q∈PmRpqRq\begin{aligned} 0 &=\mathcal{B}_p-d_pS_p+\nu_pR_p +\textstyle{\sum_{q\in\mathcal{P}}} m_{Spq}S_{q} \\ 0 &=-\left(\nu_{p}+d_{p}\right)R_{p} +\textstyle{\sum_{q\in\mathcal{P}}} m_{Rpq}R_{q} \end{aligned} 00=Bp−dpSp+νpRp+∑q∈PmSpqSq=−(νp+dp)Rp+∑q∈PmRpqRq
Want to solve for Sp,RpS_p,R_pSp,Rp. Here, it is best (crucial in fact) to remember some linear algebra. Write system in vector form:
0=b−dS+νR+MSS0=−(ν+d)R+MRR\begin{aligned} \mathbf{0} &=\mathbf{b}-\mathbf{d}\mathbf{S}+\mathbf{\nu}\mathbf{R}+\mathcal{M}^S\mathbf{S} \\ \mathbf{0} &=-\left(\mathbf{\nu}+\mathbf{d}\right)\mathbf{R}+\mathcal{M}^R\mathbf{R} \end{aligned} 00=b−dS+νR+MSS=−(ν+d)R+MRR
where S,R,b∈R∣P∣\mathbf{S},\mathbf{R},\mathbf{b}\in\mathbb{R}^{|\mathcal{P}|}S,R,b∈R∣P∣, d,ν,MS,MR\mathbf{d},\mathbf{\nu},\mathcal{M}^S,\mathcal{M}^Rd,ν,MS,MR ∣P∣×∣P∣|\mathcal{P}|\times|\mathcal{P}|∣P∣×∣P∣-matrices (d,ν\mathbf{d},\mathbf{\nu}d,ν diagonal)
Recall second equation:
0=−(ν+d)R+MRR⇔(MR−ν−d)R=0\mathbf{0} =-\left(\mathbf{\nu}+\mathbf{d}\right)\mathbf{R}+\mathcal{M}^R\mathbf{R} \Leftrightarrow (\mathcal{M}^R-\mathbf{\nu}-\mathbf{d})\mathbf{R}=\mathbf{0} 0=−(ν+d)R+MRR⇔(MR−ν−d)R=0
So unique solution R=0\mathbf{R}=\mathbf{0}R=0 if MR−ν−d\mathcal{M}^R-\mathbf{\nu}-\mathbf{d}MR−ν−d invertible. Is it?
We have been here before!
From spectrum shift, s(MR−ν−d)=−minp∈P(νp+dp)<0s(\mathcal{M}^R-\mathbf{\nu}-\mathbf{d})=-\min_{p\in\mathcal{P}}(\nu_p+d_p)<0s(MR−ν−d)=−minp∈P(νp+dp)<0
So, given L=I=0\mathbf{L}=\mathbf{I}=\mathbf{0}L=I=0, R=0\mathbf{R}=\mathbf{0}R=0 is the unique equilibrium and
limt→∞R(t)=0\lim_{t\to\infty}\mathbf{R}(t)=\mathbf{0} t→∞limR(t)=0
⟹ \implies⟹ DFE has L=I=R=0\mathbf{L}=\mathbf{I}=\mathbf{R}=\mathbf{0}L=I=R=0
DFE has L=I=R=0\mathbf{L}=\mathbf{I}=\mathbf{R}=\mathbf{0}L=I=R=0 and b−dS+MSS=0\mathbf{b}-\mathbf{d}\mathbf{S}+\mathcal{M}^S\mathbf{S}=\mathbf{0}b−dS+MSS=0, i.e.,
S=(d−MS)−1b\mathbf{S}=(\mathbf{d}-\mathcal{M}^S)^{-1}\mathbf{b} S=(d−MS)−1b
Recall: −MS-\mathcal{M}^S−MS singular M-matrix. From previous reasoning, d−MS\mathbf{d}-\mathcal{M}^Sd−MS has instability modulus shifted right by minp∈Pdp\min_{p\in\mathcal{P}}d_pminp∈Pdp. So:
Second point ⟹ (d−MS)−1>0 ⟹ (d−MS)−1b>0\implies (\mathbf{d}-\mathcal{M}^S)^{-1}>\mathbf{0}\implies (\mathbf{d}-\mathcal{M}^S)^{-1}\mathbf{b}> \mathbf{0}⟹(d−MS)−1>0⟹(d−MS)−1b>0 (would have ≫0\gg\mathbf{0}≫0 if MS\mathcal{M}^SMS irreducible)
So DFE makes sense with
(S,L,I,R)=((d−MS)−1b,0,0,0)(\mathbf{S},\mathbf{L},\mathbf{I},\mathbf{R})=\left((\mathbf{d}-\mathcal{M}^S)^{-1}\mathbf{b},\mathbf{0},\mathbf{0},\mathbf{0}\right) (S,L,I,R)=((d−MS)−1b,0,0,0)
Use next generation method with Ξ={L1,…,L∣P∣,I1,…,I∣P∣}\Xi=\{L_1,\ldots,L_{|\mathcal{P}|},I_1,\ldots,I_{|\mathcal{P}|}\}Ξ={L1,…,L∣P∣,I1,…,I∣P∣}, Ξ′=F−V\Xi'=\mathcal{F}-\mathcal{V}Ξ′=F−V
F=(Φ1,…,Φ∣P∣,0,…,0)T\mathcal{F}=\left(\Phi_1,\ldots,\Phi_{|\mathcal{P}|},0,\ldots,0\right)^T F=(Φ1,…,Φ∣P∣,0,…,0)T
V=((ε1+d1)L1−∑q∈PmL1qLq⋮(ε∣P∣+d∣P∣)L∣P∣−∑q∈PmL∣P∣qLq−ε1L1+(γ1+d1)I1−∑q∈PmI1qIq⋮−ε∣P∣L∣P∣+(γ∣P∣+d∣P∣)I∣P∣−∑q∈PmI∣P∣qIq)\mathcal{V}= \begin{pmatrix} \left( \varepsilon_{1}+d_{1}\right)L_{1} -\sum\limits_{q\in\mathcal{P}} m_{L1q}L_{q} \\ \vdots \\ \left( \varepsilon_{|\mathcal{P}|}+d_{|\mathcal{P}|}\right)L_{|\mathcal{P}|} -\sum\limits_{q\in\mathcal{P}} m_{L|\mathcal{P}|q}L_{q} \\ -\varepsilon_1L_1+(\gamma_1+d_1)I_1 -\sum\limits_{q\in\mathcal{P}} m_{I1q}I_{q} \\ \vdots \\ -\varepsilon_{|\mathcal{P}|}L_{|\mathcal{P}|} +(\gamma_{|\mathcal{P}|}+d_{|\mathcal{P}|})I_{|\mathcal{P}|} -\sum\limits_{q\in\mathcal{P}} m_{I|\mathcal{P}|q}I_{q} \end{pmatrix} V=⎝⎛(ε1+d1)L1−q∈P∑mL1qLq⋮(ε∣P∣+d∣P∣)L∣P∣−q∈P∑mL∣P∣qLq−ε1L1+(γ1+d1)I1−q∈P∑mI1qIq⋮−ε∣P∣L∣P∣+(γ∣P∣+d∣P∣)I∣P∣−q∈P∑mI∣P∣qIq⎠⎞
Differentiate w.r.t. Ξ\XiΞ:
DF=(∂Φ1∂L1⋯∂Φ1∂L∣P∣∂Φ1∂I1⋯∂Φ1∂I∣P∣⋮⋮⋮⋮∂Φ∣P∣∂L1⋯∂Φ∣P∣∂L∣P∣∂Φ∣P∣∂I1⋯∂Φ∣P∣∂I∣P∣0⋯00⋯0⋮⋮⋮⋮0⋯00⋯0)D\mathcal{F} = \begin{pmatrix} \dfrac{\partial\Phi_1}{\partial L_1} & \cdots & \dfrac{\partial\Phi_1}{\partial L_{|\mathcal{P}|}} & \dfrac{\partial\Phi_1}{\partial I_1} & \cdots & \dfrac{\partial\Phi_1}{\partial I_{|\mathcal{P}|}} \\ \vdots & & \vdots & \vdots & & \vdots \\ \dfrac{\partial\Phi_{|\mathcal{P}|}}{\partial L_1} & \cdots & \dfrac{\partial\Phi_{|\mathcal{P}|}}{\partial L_{|\mathcal{P}|}} & \dfrac{\partial\Phi_{|\mathcal{P}|}}{\partial I_1} & \cdots & \dfrac{\partial\Phi_{|\mathcal{P}|}}{\partial I_{|\mathcal{P}|}} \\ 0 & \cdots & 0 & 0 & \cdots & 0 \\ \vdots & & \vdots & \vdots & & \vdots \\ 0 & \cdots & 0 & 0 & \cdots & 0 \end{pmatrix} DF=⎝⎛∂L1∂Φ1⋮∂L1∂Φ∣P∣0⋮0⋯⋯⋯⋯∂L∣P∣∂Φ1⋮∂L∣P∣∂Φ∣P∣0⋮0∂I1∂Φ1⋮∂I1∂Φ∣P∣0⋮0⋯⋯⋯⋯∂I∣P∣∂Φ1⋮∂I∣P∣∂Φ∣P∣0⋮0⎠⎞
Note that
∂Φp∂Lk=∂Φp∂Ik=0\frac{\partial\Phi_p}{\partial L_k}=\frac{\partial\Phi_p}{\partial I_k}=0 ∂Lk∂Φp=∂Ik∂Φp=0
whenever k≠pk\neq pk=p, so
DF=(diag(∂Φ1∂L1,…,∂Φ∣P∣∂L∣P∣)diag(∂Φ1∂I1,…,∂Φ∣P∣∂I∣P∣)00)D\mathcal{F} = \begin{pmatrix} \mathsf{diag}\left( \frac{\partial\Phi_1}{\partial L_1},\ldots,\frac{\partial\Phi_{|\mathcal{P}|}}{\partial L_{|\mathcal{P}|}}\right) & \mathsf{diag}\left( \frac{\partial\Phi_1}{\partial I_1},\ldots,\frac{\partial\Phi_{|\mathcal{P}|}}{\partial I_{|\mathcal{P}|}}\right) \\ \mathbf{0} & \mathbf{0} \end{pmatrix} DF=(diag(∂L1∂Φ1,…,∂L∣P∣∂Φ∣P∣)0diag(∂I1∂Φ1,…,∂I∣P∣∂Φ∣P∣)0)
If Φp=βpSpIp\Phi_p=\beta_pS_pI_pΦp=βpSpIp, then
If Φp=βpSpIpNp\Phi_p=\beta_p\dfrac{S_pI_p}{N_p}Φp=βpNpSpIp, then
In both cases, ∂/∂L\partial/\partial L∂/∂L block is zero so
F=DF(DFE)=(0diag(∂Φ1∂I1,…,∂Φ∣P∣∂I∣P∣)00)F=D\mathcal{F}(DFE)= \begin{pmatrix} \mathbf{0} & \mathsf{diag}\left( \frac{\partial\Phi_1}{\partial I_1},\ldots,\frac{\partial\Phi_{|\mathcal{P}|}}{\partial I_{|\mathcal{P}|}}\right) \\ \mathbf{0} & \mathbf{0} \end{pmatrix} F=DF(DFE)=(00diag(∂I1∂Φ1,…,∂I∣P∣∂Φ∣P∣)0)
V=(diagp(εp+dp)−ML0−diagp(εp)diagp(γp+dp)−MI)V= \begin{pmatrix} \mathsf{diag}_p(\varepsilon_p+d_p)-\mathcal{M}^L & \mathbf{0} \\ -\mathsf{diag}_p(\varepsilon_p) & \mathsf{diag}_p(\gamma_p+d_p)-\mathcal{M}^I \end{pmatrix} V=(diagp(εp+dp)−ML−diagp(εp)0diagp(γp+dp)−MI)
where diagp(zp)=diag(z1,…,z∣P∣)\mathsf{diag}_p(z_p)=\mathsf{diag}(z_1,\ldots,z_{|\mathcal{P}|})diagp(zp)=diag(z1,…,z∣P∣). Inverse of VVV easy (2×22\times 22×2 block lower triangular):
V−1=((diagp(εp+dp)−ML)−10V~21−1(diagp(γp+dp)−MI)−1)V^{-1} = \begin{pmatrix} \left(\mathsf{diag}_p(\varepsilon_p+d_p)-\mathcal{M}^L\right)^{-1} & \mathbf{0} \\ \tilde V_{21}^{-1} & \left(\mathsf{diag}_p(\gamma_p+d_p)-\mathcal{M}^I\right)^{-1} \end{pmatrix} V−1=((diagp(εp+dp)−ML)−1V~21−10(diagp(γp+dp)−MI)−1)
where
V~21−1=(diagp(εp+dp)−ML)−1diagp(εp)(diagp(γp+dp)−MI)−1\tilde V_{21}^{-1}= \left(\mathsf{diag}_p(\varepsilon_p+d_p)-\mathcal{M}^L\right)^{-1} \mathsf{diag}_p(\varepsilon_p) \left(\mathsf{diag}_p(\gamma_p+d_p)-\mathcal{M}^I\right)^{-1} V~21−1=(diagp(εp+dp)−ML)−1diagp(εp)(diagp(γp+dp)−MI)−1
Next generation matrix
FV−1=(0F1200)(V~11−10V~21−1V~22−1)=(F12V~21−1F12V~22−100)FV^{-1}= \begin{pmatrix} \mathbf{0} & F_{12} \\ \mathbf{0} & \mathbf{0} \end{pmatrix} \begin{pmatrix} \tilde V_{11}^{-1} & \mathbf{0} \\ \tilde V_{21}^{-1} & \tilde V_{22}^{-1} \end{pmatrix} = \begin{pmatrix} F_{12}\tilde V_{21}^{-1} & F_{12}\tilde V_{22}^{-1} \\ \mathbf{0} & \mathbf{0} \end{pmatrix} FV−1=(00F120)(V~11−1V~21−10V~22−1)=(F12V~21−10F12V~22−10)
where V~ij−1\tilde V_{ij}^{-1}V~ij−1 is block ijijij in V−1V^{-1}V−1. So
R0=ρ(F12V~21−1)\mathcal{R}_0=\rho\left(F_{12}\tilde{V}_{21}^{-1}\right) R0=ρ(F12V~21−1)
i.e.,
R0=ρ(diag(∂Φ1∂I1,…,∂Φ∣P∣∂I∣P∣)(diagp(εp+dp)−ML)−1diagp(εp)(diagp(γp+dp)−MI)−1)\mathcal{R}_0=\rho\Biggl( \mathsf{diag}\left( \frac{\partial\Phi_1}{\partial I_1},\ldots,\frac{\partial\Phi_{|\mathcal{P}|}}{\partial I_{|\mathcal{P}|}}\right) \left(\mathsf{diag}_p(\varepsilon_p+d_p)-\mathcal{M}^L\right)^{-1} \\ \qquad\qquad\qquad\qquad\qquad\qquad \mathsf{diag}_p(\varepsilon_p) \left(\mathsf{diag}_p(\gamma_p+d_p)-\mathcal{M}^I\right)^{-1} \Biggr) R0=ρ(diag(∂I1∂Φ1,…,∂I∣P∣∂Φ∣P∣)(diagp(εp+dp)−ML)−1diagp(εp)(diagp(γp+dp)−MI)−1)
Define R0\mathcal{R}_0R0 for the ∣P∣|\mathcal{P}|∣P∣-SLIRS as
Then the DFE
is locally asymptotically stable if R0<1\mathcal{R}_0<1R0<1 and unstable if R0>1\mathcal{R}_0>1R0>1
From PvdD & Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Bulletin of Mathematical Biology 180(1-2): 29-48 (2002)
JA & S Portet. Epidemiological implications of mobility between a large urban centre and smaller satellite cities. Journal of Mathematical Biology 71(5):1243-1265 (2015)
Winnipeg as urban centre and 3 smaller satellite cities: Portage la Prairie, Selkirk and Steinbach
Assume myxm_{yx}myx movement rate from city xxx to city yyy. Ceteris paribus, Nx′=−myxNxN_x'=-m_{yx}N_xNx′=−myxNx, so Nx(t)=Nx(0)e−myxtN_x(t)=N_x(0)e^{-m_{yx}t}Nx(t)=Nx(0)e−myxt. Therefore, after one day, Nx(1)=Nx(0)e−myxN_x(1)=N_x(0)e^{-m_{yx}}Nx(1)=Nx(0)e−myx, i.e.,
myx=−ln(Nx(1)Nx(0))m_{yx}=-\ln\left(\frac{N_x(1)}{N_x(0)}\right) myx=−ln(Nx(0)Nx(1))
Now, Nx(1)=Nx(0)−TyxN_x(1)=N_x(0)-T_{yx}Nx(1)=Nx(0)−Tyx, where TyxT_{yx}Tyx number of individuals going from xxx to yyy / day. So
myx=−ln(1−TyxNx(0))m_{yx}=-\ln\left(1-\frac{T_{yx}}{N_x(0)}\right) myx=−ln(1−Nx(0)Tyx)
Computed for all pairs (W,i)(W,i)(W,i) and (i,W)(i,W)(i,W) of cities
with disease: R0x=1.5\mathcal{R}_0^x=1.5R0x=1.5; without disease: R0x=0.5\mathcal{R}_0^x=0.5R0x=0.5. Each box and corresponding whiskers are 10,000 simulations
PLP and Steinbach have comparable populations but with parameters used, only PLP can cause the general R0\mathcal{R}_0R0 to take values larger than 1 when R0W<1\mathcal{R}_0^W<1R0W<1
This is due to the movement rate: if M=0\mathcal{M}=0M=0, then
R0=max{R0W,R01,R02,R03},\mathcal{R}_0=\max\{\mathcal{R}_0^W,\mathcal{R}_0^1,\mathcal{R}_0^2,\mathcal{R}_0^3\}, R0=max{R0W,R01,R02,R03},
since FV−1FV^{-1}FV−1 is then block diagonal
Movement rates to and from PLP are lower →\rightarrow→ situation closer to uncoupled case and R01\mathcal{R}_0^1R01 has more impact on the general R0\mathcal{R}_0R0
Plots as functions of R01\mathcal{R}_0^1R01 in PLP and the reduction of movement between Winnipeg and PLP. Left: general R0\mathcal{R}_0R0. Right: Attack rate in Winnipeg
Given
skp′=fkp(Sp,Ip)iℓp′=gℓp(Sp,Ip)\begin{aligned} s_{kp}' &= f_{kp}(S_p,I_p) \\ i_{\ell p}' &= g_{\ell p}(S_p,I_p) \end{aligned} skp′iℓp′=fkp(Sp,Ip)=gℓp(Sp,Ip)
with known properties, what is known of
skp′=fkp(Sp,Ip)+∑q∈Pmkpqskqiℓp′=gℓp(Sp,Ip)+∑q∈Pmℓpqiℓq?\begin{aligned} s_{kp}' &= f_{kp}(S_p,I_p)+\textstyle{\sum_{q\in\mathcal{P}}} m_{kpq}s_{kq} \\ i_{\ell p}' &= g_{\ell p}(S_p,I_p)+\textstyle{\sum_{q\in\mathcal{P}}} m_{\ell pq}i_{\ell q} \qquad ? \end{aligned} skp′iℓp′=fkp(Sp,Ip)+∑q∈Pmkpqskq=gℓp(Sp,Ip)+∑q∈Pmℓpqiℓq?
YES, coupling together backward bifurcating units can lead to a system-level backward bifurcation
JA, Ducrot & Zongo. A metapopulation model for malaria with transmission-blocking partial immunity in hosts. Journal of Mathematical Biology 64(3):423-448 (2012)
"Converse" problem to inheritance problem. Given
with known properties, does
skp′=fkp(Sp,Ip)+∑q∈Pmkpqskqiℓp′=gℓp(Sp,Ip)+∑q∈Pmℓpqiℓq\begin{aligned} s_{kp}' &= f_{kp}(S_p,I_p)+\textstyle{\sum_{q\in\mathcal{P}}} m_{kpq}s_{kq} \\ i_{\ell p}' &= g_{\ell p}(S_p,I_p)+\textstyle{\sum_{q\in\mathcal{P}}} m_{\ell pq}i_{\ell q} \end{aligned} skp′iℓp′=fkp(Sp,Ip)+∑q∈Pmkpqskq=gℓp(Sp,Ip)+∑q∈Pmℓpqiℓq
exhibit some behaviours not observed in the uncoupled system?
E.g.: units have {R0<1 ⟹ \{\mathcal{R}_0<1\implies{R0<1⟹ DFE GAS, R0>1 ⟹ \mathcal{R}_0>1\impliesR0>1⟹ 1 GAS EEP}\}} behaviour, metapopulation has periodic solutions
Can there be situations where some patches are at the DFE and others at an EEP?
This is the problem of mixed equilibria
This is a metapopulation-specific problem, not one of inheritance of dynamical properties!
[Patch level] Patch p∈Pp\in\mathcal{P}p∈P at equilibrium is empty if Xp⋆=0X_p^\star=0Xp⋆=0, at the disease-free equilibrium if Xp⋆=(sk1p⋆,…,skup⋆,0,…,0)X_p^\star=(s_{k_1p}^\star,\ldots,s_{k_up}^\star,0,\ldots,0)Xp⋆=(sk1p⋆,…,skup⋆,0,…,0), where k1,…,kuk_1,\ldots,k_uk1,…,ku are some indices with 1≤u≤∣U∣1\leq u\leq|\mathcal{U}|1≤u≤∣U∣ and sk1p⋆,…,skup⋆s_{k_1p}^\star,\ldots,s_{k_up}^\starsk1p⋆,…,skup⋆ are positive, and at an endemic equilibrium if Xp≫0X_p\gg 0Xp≫0
[Metapopulation level] A population-free equilibrium has all patches empty. A metapopulation disease-free equilibrium has all patches at the disease-free equilibrium for the same compartments. A metapopulation endemic equilibrium has all patches at an endemic equilibrium
A mixed equilibrium is an equilibrium such that
E.g.,
((S1,I1,R1),(S2,I2,R2))=((+,0,0),(+,+,+))((S_1,I_1,R_1),(S_2,I_2,R_2))=((+,0,0),(+,+,+)) ((S1,I1,R1),(S2,I2,R2))=((+,0,0),(+,+,+))
is mixed, so is
((S1,I1,R1),(S2,I2,R2))=((+,0,0),(+,0,+))((S_1,I_1,R_1),(S_2,I_2,R_2))=((+,0,0),(+,0,+)) ((S1,I1,R1),(S2,I2,R2))=((+,0,0),(+,0,+))
Suppose that movement is similar for all compartments (MSAC) and that the system is at equilibrium
Note that MSAC ⟹ \implies⟹ Ac=A\mathcal{A}^c=\mathcal{A}Ac=A and Dc=D\mathcal{D}^c=\mathcal{D}Dc=D for all c∈Cc\in\mathcal{C}c∈C
More is needed on inheritance problem, in particular GAS part (Li, Shuai, Kamgang, Sallet, and older stuff: Michel & Miller, Šiljak)
Incorporate travel time (delay) and events (infection, recovery, death ..) during travel
What is the minimum complexity of the movement functions mmm below
skp′=fkp(Sp,Ip)+∑q∈Pmkpq(S,I)skqiℓp′=gℓp(Sp,Ip)+∑q∈Pmℓpq(S,I)iℓq\begin{aligned} s_{kp}' &= f_{kp}(S_p,I_p)+\textstyle{\sum_{q\in\mathcal{P}}} m_{kpq}(S,I)s_{kq} \\ i_{\ell p}' &= g_{\ell p}(S_p,I_p)+\textstyle{\sum_{q\in\mathcal{P}}} m_{\ell pq}(S,I)i_{\ell q} \end{aligned} skp′iℓp′=fkp(Sp,Ip)+∑q∈Pmkpq(S,I)skq=gℓp(Sp,Ip)+∑q∈Pmℓpq(S,I)iℓq
required to observe a metapopulation-induced behaviour?
SLIAR_metapop_rhs <- function(t, x, p) { with(as.list(x), { S = x[p$idx_S] L = x[p$idx_L] I = x[p$idx_I] A = x[p$idx_A] R = x[p$idx_R] N = S + L + I + A + R Phi = p$beta * S * (I + p$eta * A) dS = S - Phi + p$M %*% S dL = Phi - p$epsilon * L + p$M %*% L dI = (1 - p$pi) * p$epsilon * L - p$gammaI * I + p$M %*% I dA = p$pi * p$epsilon * L - p$gammaA * A + p$M %*% A dR = p$gammaI * I + p$gammaA * A + p$M %*% R list(c(dS, dL, dI, dA, dR)) }) }
pop = c(34.017, 1348.932, 1224.614, 173.593, 93.261) * 1e+06 countries = c("Canada", "China", "India", "Pakistan", "Philippines") T = matrix(data = c(0, 1268, 900, 489, 200, 1274, 0, 678, 859, 150, 985, 703, 0, 148, 58, 515, 893, 144, 0, 9, 209, 174, 90, 2, 0), nrow = 5, ncol = 5, byrow = TRUE)
p = list() # Use the approximation explained in Arino & Portet (JMB 2015) p$M = mat.or.vec(nr = dim(T)[1], nc = dim(T)[2]) for (from in 1:5) { for (to in 1:5) { p$M[to, from] = -log(1 - T[from, to]/pop[from]) } p$M[from, from] = 0 } p$M = p$M - diag(colSums(p$M))
p$P = dim(p$M)[1] p$idx_S = 1:p$P p$idx_L = (p$P + 1):(2 * p$P) p$idx_I = (2 * p$P + 1):(3 * p$P) p$idx_A = (3 * p$P + 1):(4 * p$P) p$idx_R = (4 * p$P + 1):(5 * p$P) p$eta = rep(0.3, p$P) p$epsilon = rep((1/1.5), p$P) p$pi = rep(0.7, p$P) p$gammaI = rep((1/5), p$P) p$gammaA = rep((1/3), p$P) # The desired values for R_0. Here we take something simple R_0 = rep(1.5, p$P)
# Set initial conditions. For example, we start with 2 # infectious individuals in Canada. L0 = mat.or.vec(p$P, 1) I0 = mat.or.vec(p$P, 1) A0 = mat.or.vec(p$P, 1) R0 = mat.or.vec(p$P, 1) I0[1] = 2 S0 = pop - (L0 + I0 + A0 + R0) # Vector of initial conditions to be passed to ODE solver. IC = c(S = S0, L = L0, I = I0, A = A0, R = R0) # Time span of the simulation (5 years here) tspan = seq(from = 0, to = 5 * 365.25, by = 0.1)
Let us take R0=1.5\mathcal{R}_0=1.5R0=1.5 for patches in isolation. Solve R0\mathcal{R}_0R0 for β\betaβ
β=R0S(0)(1−πpγIp+πpηpγAp)−1\beta=\frac{\mathcal{R}_0}{S(0)} \left( \frac{1-\pi_p}{\gamma_{Ip}} +\frac{\pi_p\eta_p}{\gamma_{Ap}} \right)^{-1} β=S(0)R0(γIp1−πp+γApπpηp)−1
for (i in 1:p$P) { p$beta[i] = R_0[i] / S0[i] * 1/((1 - p$pi[i])/p$gammaI[i] + p$pi[i] * p$eta[i]/p$gammaA[i]) }
# Call the ODE solver sol <- deSolve::ode(y = IC, times = tspan, func = SLIAR_metapop_rhs, parms = p) ## DLSODA- At current T (=R1), MXSTEP (=I1) steps ## taken on this call before reaching TOUT ## In above message, I1 = 5000 ## ## In above message, R1 = 117.498
The output I copy above means the integration went wrong. The problem is the sie difference between countries, in particular China and Canada..
Need to play with movement rates and initial conditions. Will not explain here
Simulating the ODE metapopulation as a continuous time Markov chain (CTMC) is a little more complicated
However, still relatively easy to do if being careful using R libraries such as adaptivetau or GillespieSSA2
R
adaptivetau
GillespieSSA2
Space is a fundamental component of the epidemic spread process and cannot be ignored, both in modelling and in public health decision making
One way to model space is to use metapopulation models
Metapopulation models are easy to analyse locally, give interesting problems at the global level and are easy to simulate
Simulation (deterministic and stochastic) can be costly in RAM and cycles
Metapopulation models are not the only solution
<div style = "text-align: justify">
</div>