Supporting Information for
Spotting epidemic keystones by 𝑹
𝟎sensitivity analysis: High-risk stations in the Tokyo metropolitan area
Kenta Yashima1, 2, Akira Sasaki1, 3
1 Department of Evolutionary Studies of Biosystems, the Graduate University for Advanced Studies (SOKENDAI), Hayama, Kanagawa, Japan.
2Meiji Institute for Advanced Study of Mathematical Sciences, Meiji University, Nakano, Tokyo, Japan.
3Evolution and Ecology Program, International Institute for Applied Systems Analysis, Laxenburg, Austria.
Supporting Information S1. Derivation of the next generation matrix 𝑳
Derivation of the next generation matrix 𝑳 for an infectious disease in a metropolitan area is given. By applying a linear approximation 𝑆!! 𝑡 ≅𝑁!! and 𝑆!"! 𝑡 ≅𝑁!"! to the integral forms of Eqs. (2), (5), we have the following renewal equations.
𝐼!! 𝑡 =𝛽 !𝑒!!" 2𝐼!! 𝑡−𝜏 +𝐼!! 𝑡−𝜏 𝑁!!𝑑𝜏
!
(A1)
𝐼!"! 𝑡 = 𝛽 !𝑒!!" 𝐼!! 𝑡−𝜏 +𝐼!! 𝑡−𝜏
!
+𝐼!! 𝑡−𝜏 𝑁!"!𝑑𝜏
(A2)
Here, 𝐼!! 𝑡 ≡ !𝐼!"! 𝑡 and 𝐼!! 𝑡 ≡ !𝐼!"! 𝑡 denote the number of infectious commuters in i-th home population and in j-th work population, respectively. By summing over indices i and j of Eq. (A2) and introducing a vector notation for non-commuting resident population (𝑰𝑹 𝑡 ≡ 𝐼!! 𝑡 ,𝐼!! 𝑡 ,⋯𝐼!! 𝑡 !), commuting home population (𝑰𝑯 𝑡 ≡ 𝐼!! 𝑡 ,𝐼!! 𝑡 ,⋯𝐼!! 𝑡 !), and commuting work population (𝑰𝑾 𝑡 ≡ 𝐼!! 𝑡 ,𝐼!! 𝑡 ,⋯𝐼!! 𝑡 !), Eqs. (A1), (A2) can be summarized in the 3 × 3 block matrix form as
𝑰𝑹 𝑡 𝑰𝑯 𝑡 𝑰𝑾 𝑡
=𝛽 𝑒!!"
𝑻!! 𝑻!" 0
𝑻!" 𝑻!! 𝑻!"
𝑻!" 𝑻!" 𝑻!!
𝑰𝑹 𝑡−𝜏 𝑰𝑯 𝑡−𝜏 𝑰𝑾 𝑡−𝜏
! 𝑑𝜏
! . (A3)
Here, element of this block matrix 𝑻!" is a 𝑀×𝑀 matrix and denotes the transmission from type n population to type m population (𝑚,𝑛 ∈ 𝑅,𝐻,𝑊 , R:
non-commuting resident population, H: commuting home population, W: commuting work population), where each element of the matrices are given as 𝑻!! !" ≡ 2𝑁!!𝛿!", 𝑻!" !" ≡ 𝑁!!𝛿!", 𝑻!" !" ≡𝑁!!𝛿!", 𝑻!! !" ≡ 𝑁!!𝛿!", 𝑻!" !" ≡𝑁!", 𝑻!" !" ≡
𝑁!", 𝑻!" !" ≡𝑁!", and 𝑻!! !" ≡ 𝑁!!𝛿!" (𝟎: 𝑀×𝑀 zero matrix). Here, it should
be noted that the matrix 𝑻!" includes the information about the host population structure only and the epidemiological information is not included. Since 𝛽 !!𝑒!!"𝑻!"𝑰𝒏 𝑡−𝜏 𝑑𝜏 gives the number of infectious individual from type n population to type m population, the asymptotic ratio between these population types at exponential growth phase are given as 𝛽 !!𝑒!!"𝑻!!𝑑𝜏= 𝛽/𝛾 𝑻!". Given this the next generation matrix 𝑳 is given as a 3×3 block matrix form as
𝑳 =𝛽 𝑒!!"
𝑻!! 𝑻!" 0
𝑻!" 𝑻!! 𝑻!"
𝑻!" 𝑻!" 𝑻!! 𝑑𝜏
!
!
= 𝛽 𝛾
𝑻!! 𝑻!" 0
𝑻!" 𝑻!! 𝑻!"
𝑻!" 𝑻!" 𝑻!! .
(A4)
Here, each element of the block matrix gives the asymptotic ratio between different population types. Accordingly, the dominant eigenvalue of this next generation matrix 𝑳 gives the basic reproductive ratio 𝑅! =𝜌 𝑳 (𝜌 ⋯ : spectral radius).
For the dominant eigenvalue 𝑅!, the elements of corresponding left and right eigenvector give the reproductive value and the relative ratio of exponentially growing infected population, respectively [1–3]. The reproductive value (i.e., element of left eigenvector 𝑣!!, 𝑣!!,𝑣!! ) of each local population is given in Figure S2-A as a function of its local population size 𝑁!!, 𝑁!!,𝑁!! . The relative ratio of exponentially growing infected population (i.e., element of right eigenvector 𝑤!!, 𝑤!!,𝑤!! ) of each local population is given in Figure S2-B as a function of its local population size
𝑁!!, 𝑁!!,𝑁!! . Both of the values increase as the local population size increase. This means that, a local population with a larger population size has a larger impact on the overall epidemic dynamics and also has a higher risk of infection. Furthermore, the most notable point is that, the results can be clustered into two distinct groups. This separation can be explained in the relation to the largest work population (i.e., working/studying area of Shinjuku station). For the results of non-commuting resident population and commuting home population (Figure S2-A2, 3, S2-B2, 3), the local populations in the upper cluster has at least one commuter who is working/studying at the largest work population. The horizontally layered colored structure can be clearly explained by the number of commuters to the largest work population. On the other hand, no one from the local populations in the lower cluster is working/studying at the largest work population. For the results of commuting population at work population (Figure S2-A1, S2-B1), the upper cluster is the largest work population itself and the lower cluster is consisted from other local populations. The clear distinction observed in the change in the basic reproductive ratio 𝛿𝑅! (see Figure 3), can be attributed to this distinction in the eigenvectors.
Supporting Information S2. Final size of epidemic
Overall damage can be evaluated by the global final size of epidemic 𝛹, which is defined as a ratio of infected individuals whom has ever acquired infection during the epidemic period. For this we define the local final size of epidemic 𝛹!! within the non-commuting resident population at i-th station as
𝛹!! ≡𝑅!! ∞
𝑁!! = 𝑁!!−𝑆!! ∞
𝑁!! ≅1−𝑆!! ∞
𝑆!! 0 , (A5)
and the local final size of epidemic 𝛹!"! within the commuting population consisted from a commuters residing at i-th home population and working at j-th work population as
𝛹!"! ≡ 𝑅!"! ∞
𝑁!"! =𝑁!"!−𝑆!"! ∞
𝑁!"! ≅1−𝑆!"! ∞
𝑆!"! 0 . (A6)
Here, we have used the approximations 𝑆!! 0 ≅ 𝑁!! and 𝑆!"! 0 ≅𝑁!"!. These can be calculated in the following way. The integral forms of Eqs. (1), (4) and Eqs. (3), (6) are given as following.
𝑆!! 𝑡 = 𝑆!! 0 exp −𝛽 2𝐼!! 𝑠 + 𝐼!"!
!
𝑠 𝑑𝑠
!
!
(A7)
𝑆!"! 𝑡 =𝑆!"! 0 exp −𝛽 𝐼!! 𝑠
!
!
+ 𝐼!"! 𝑠 +𝐼!"! 𝑠
!
𝑑𝑠
(A8)
𝑅!! 𝑡 =𝛾 !𝐼!! 𝑡 𝑑𝑠
!
(A9)
𝑅!"! 𝑡 =𝛾 !𝐼!"! 𝑡 𝑑𝑠
!
(A10)
Substitution of Eqs. (A7)-(A10) with 𝑡= ∞ to the definition of the local final size of epidemic Eq. (A5) and Eq. (A6) yields
𝛹!! =1− exp −𝛽 2𝐼!! 𝑠 + 𝐼!"!
!
𝑠 𝑑𝑠
!
!
=1− exp −𝛽
𝛾 2𝑁!!𝛹!! + 𝑁!"!𝛹!"!
!
,
(A11)
𝛹!"! =1− exp −𝛽 ! 𝐼!! 𝑠
!
+ 𝐼!"! 𝑠 +𝐼!"! 𝑠
!
𝑑𝑠
=1− exp −𝛽
𝛾 𝑁!!𝛹!! + 𝑁!"!𝛹!"! +𝑁!"!𝛹!"!
!
.
(A12)
This system of transcendental equation with 𝑀 𝑀+1 ~100,000 equations can be solved numerically; by recursive calculation starting from 𝛹!! = 1 and 𝛹!"! = 1 until all the values converges to a fixed point. Once this value 𝛹!"! has been obtained, the local final size of epidemic for commuting population at each home population 𝛹!! and each work population 𝛹!! can be obtained in the following way.
𝛹!! ≡ 𝑅!! ∞
𝑁!! = !𝑁!"!𝛹!"!
𝑁!! (A13)
𝛹!! ≡ 𝑅!! ∞
𝑁!! = !𝑁!"!𝛹!"!
𝑁!! (A14)
Here the number of recovered individuals at i-th home population (j-th work population) is denoted as 𝑅!! 𝑡 ≡ !𝑅!"! 𝑡 (𝑅!! 𝑡 ≡ !𝑅!"! 𝑡 ). The final size of epidemic within a total non-commuting population 𝛹! and total commuting population 𝛹! can be obtained in similar fashion from 𝛹!! and 𝛹!"! as
𝛹! = !𝑅!! ∞ 𝑁!!
! = !𝑁!!𝛹!!
𝑁! , (A15)
𝛹! = ! !𝑅!"! ∞
𝑁!"!
!
!
= ! !𝑁!"!𝛹𝑖𝑗!
𝑁! . (A16)
Then the global final size of epidemic 𝛹 is given from 𝛹! and 𝛹! as
𝛹≡ !𝑅!! ∞ + ! !𝑅!"! ∞
𝑁 =𝑁!𝛹! +𝑁!𝛹!
𝑁 . (A17)
Dependence of the final size of epidemic within the non-commuting population 𝛹!, the commuting population 𝛹! and for the total population 𝛹 on the infection rate 𝛽 is given in Figure S1-A. For both cases of 𝑟= 0 and 𝑟= 1, when the infection rate is small, the final sizes of epidemic are negligible such that the initial extinction of disease occurs, however as the infection rate increases the global final sizes of epidemic monotonically increase until it saturates to one. The infection rates at these disease invasion thresholds agree with 𝛽! obtained from the calculation of basic reproductive ratio (cf. Eq. 8). These critical infection rate 𝛽! was slightly larger for 𝑟= 0 (𝛽! =9.210485×10!!) compare to 𝑟=1 (𝛽! = 9.207523×10!!), this can be attributed to the larger total population size for 𝑟=1. This fact suggests that the effect of non-commuting population is minimal for the disease invasion condition, such that even thought the total population is doubled the effect to the threshold value 𝛽! is minimal. For 𝑟= 1, the final size of epidemic is larger for commuting population (𝛹!) compare to that of the non-commuting population (𝛹!), and the disease invasion occurs from slightly smaller infection rate. These can be ascribed to the fact that the commuting individuals have a higher risk of encountering infectious individuals at the
work population compare to the non-commuting residents. In such way, there is a quantitative difference between the results of 𝑟= 0 and 𝑟=1. However, qualitative differences are minor, therefore, throughout this study we have used 𝑟= 1 in the analysis.
References for supplementary materials:
1. Diekmann O, Heesterbeek JAP. Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation (Wiley Series in Mathematical & Computational Biology). JOHN WILEY & SON; 2000.
2. Caswell H. Matrix Population Models. Second. Sinauer Associates Inc.; 2001.
3. Ellner SP, Guckenheimer J. Dynamic Models in Biology. Princeton University Press; 2011.
4. Klemm K, Serrano MA, Eguiluz VM, Miguel MS. A measure of individual role in collective dynamics: spreading at criticality. Sci Rep. 2012;2: 292.
doi:10.1038/srep00292
Supporting Information S3. Supporting Figures
Figure S1: The final size of epidemic for infectious disease spread in the Tokyo metropolitan area. (A) The global final size 𝛹 of the epidemic for 𝑟=0 (blue line) and 𝑟=1 (red lines) are plotted against the infection rate 𝛽 (𝑟: ratio of non-commuting individuals to commuting individuals, see Methods section for details).
For 𝑟= 0, where all the population would commute, the result for the commuting population 𝛹! (in this case the same as the result for the total population) is only present as a blue solid line. For 𝑟=1, the type of line indicates the result for the commuting population 𝛹! (red dotted line), non-commuting population 𝛹! (red dashed line), and total population 𝛹 (red solid line), respectively. Note that the result for the total population with 𝑟= 1 (red solid line) overlaps with the result for commuting population with 𝑟=0 (blue solid lines) and is not visible on the figure. As the infection rate 𝛽 exceeds a threshold value, the global final size Ψ of the epidemic becomes non-zero and increases along with the infection rate, for both values of 𝑟.
According to the analysis of the basic reproductive ratio, the threshold value of infection 𝛽! is given as 𝛽! =9.210485×10!! for 𝑟=0 and 𝛽! =9.207523×
10!! for 𝑟= 1. These values are in good agreement with the results obtained for the final size of the epidemic. (B) The local final size of the epidemic at work population (𝛹!!), home population (𝛹!!), and resident population (𝛹!!) of each station are plotted against its local population size in figures B1, 2, and 3, respectively. The results for different infection rates are denoted by different colors, here 𝑟= 1 is used for the calculation. There is a sigmoidal dependence of local final size of epidemic on its
(A) Global final size (B1) Local final size (Work) (B2) Local final size (Home) (B3) Local final size (Resident)
population size, such that the local final size is small when the population size is small and as the population size becomes larger it will increase until it saturates to one at the larger limit. Here the location of the steep transition point will shift to the smaller side as the infection rate becomes larger. This point will become relevant in relation to the effect of countermeasures on the final size of epidemic (see Figure 4B, C).
Figure S2: The left and right eigenvectors corresponding to the dominant eigenvalue of next generation matrix for infectious disease spread in the Tokyo metropolitan area. (A) The element of the left eigenvector (𝑣!!, 𝑣!!, 𝑣!!) that gives the reproductive value of infection at each local population is plotted against its local population size (𝑁!!, 𝑁!!, 𝑁!!), where each dot represents a single station. The
“dynamic influence” introduced by Klemm et al. [4] corresponds to this value, except that they have calculated the eigenvector of the Jacobian matrix and not the next generation matrix. (B) The element of right eigenvector (𝑤!!, 𝑤!!, 𝑤!!) that gives the relative fraction of infected individuals at each local population in an exponentially growing phase is plotted against its local population size (𝑁!!, 𝑁!!, 𝑁!!), where each dot represents a single station. Results for the commuting population at each work population and home population are given in (A1, B1) and (A2, B2), respectively and the results for the non-commuting resident population at each station are given in (A3, B3). The color of each dot shows their relationship with the largest work population (Shinjuku station). Black diamonds marked with a red circle in (A1) and (B1) correspond to the largest work population, and other work populations are represented by black triangles. Each colored dot in (A2, 3) and (B2, 3) corresponds to a station that
(A2) Left eigenvector (Home)
(B2) Right eigenvector (Home)
(A3) Left eigenvector (Resident)
(A1) Left eigenvector (Work)
(B3) Right eigenvector (Resident)
(B1) Right eigenvector (Work)
has at least one commuter that travels to the largest work population and the color indicates the number of commuters who go there. Black dots correspond to stations with no commuters to the largest work population. For both commuting and non-commuting populations, the elements of the leading left and right eigenvectors were separated into two distinct groups, which can be interpreted from their relationship with the largest work population. The strong dependence of the 𝑅!-centrality on the Shinjuku station originates from this characteristic.
Figure S3: The 𝑹𝟎-centrality for every commuting pathway and residential station in the Tokyo metropolitan area (the same data as Figure 2 presented in relation to the second largest work population). The 𝑅!-centrality for each commuting population (each dot in Figure S3A corresponds to a single commuting pathway) and non-commuting population (each dot in Figure S3B corresponds to a single residential station) are given in accordance with the relation to the working population at Tokyo station. The schematic illustration above each panel describes its relationship. The 𝑅!-centralities in the commuting populations (A-1), those who commute directly to Tokyo station, (A-2), those who do not commute to Tokyo station but share a common resident station with them, (A-3): neither of them, are plotted against the population size of its working population (𝑁!!). Similarly, the 𝑅!-centralities of non-commuting population (B-1), those residing at the station area from which at least one commutes to Tokyo station, and (B-2), those residing at the station area from which no one commutes to Tokyo station, are plotted against the population size of its resident population (𝑁!!). The color of dots indicates the number of commuters to the working population at Tokyo station.
2 Home
Work Resident
2 Home
Work Resident
2 Home
Work Resident
2 Home
Work Resident
(A) Commuting population (B) Non-commuting population
(A-1) (A-2) (A-3) (B-1)
Commute to Tokyo (2nd largest)
Commute to Shinjuku (largest)
Figure S4: The 𝑹𝟎-centrality for every commuting pathway and residential station in the Tokyo metropolitan area, after vaccinating/quarantining every individual from the largest working population at Shinjuku station. The 𝑅!-centrality for each commuting population and non-commuting population after the vaccinating/quarantining all the individual from the largest working population at Shinjuku station, are given in accordance with the relation to the working population at
Commute to Tokyo (largest) 2 Home
Work Resident
2 Home
Work Resident
2 Home
Work Resident
2 Home
Work Resident
2 Home
Work Resident
(A) Commuting population (B) Non-commuting population
(A-1) (A-2) (A-3) (B-1) (B-2)
Tokyo
3 Home
Work Resident
3 Home
Work Resident
3 Home
Work Resident
3 Home
Work Resident
3 Home
Work Resident
(C) Commuting population (D) Non-commuting population
(C-1) (C-2) (C-3) (D-1) (D-2)
Commute to Shibuya (2nd largest)
Shibuya
Commute to Tokyo (largest)
Tokyo station (currently the largest susceptible work population after the removal of Shinjuku station) and Shibuya station (currently the second largest susceptible work population) are given in (A for commuting population, B for non-commuting population) and (C for commuting population, D for non-commuting population), respectively. The schematic illustration above each panel describes its relationship. The color of dots indicates the number of susceptible commuters to the working population at Tokyo station in (A, B) and to the working population at Shibuya station in (C, D).
Figure S5: The 𝑹𝟎-centrality for every commuting pathway and residential station in the Tokyo metropolitan area, after vaccinating/quarantining 1,700 individuals from the largest working population at Shinjuku station. The 𝑅!-centrality for each commuting population and non-commuting population after vaccinating/quarantining 1,700 individuals from the largest working population at Shinjuku station, are given in accordance with the relation to the working population at Shinjuku station (the largest
(A) Commuting population (B) Non-commuting population 1 Home
Work Resident
1 Home
Work Resident
1 Home
Work Resident
1 Home
Work Resident
1 Home
Work Resident
(A-1) (A-2) (A-3) (B-1) (B-2)
Shinjuku
Commute to Shinjuku (largest)
Commute to Tokyo (2nd largest)
2 Home
Work Resident
2 Home
Work Resident
2 Home
Work Resident
2 Home
Work Resident
2 Home
Work Resident
(C) Commuting population (D) Non-commuting population
(C-1) (C-2) (C-3) (D-1) (D-2)
Tokyo
Commute to Shinjuku (largest)
Commute to Tokyo (2nd largest)
susceptible work population) and Tokyo station (the second largest susceptible work population) are given in (A for commuting population, B for non-commuting population) and (C for commuting population, D for non-commuting population), respectively. The schematic illustration above each panel describes its relationship. The color of dots indicates the number of susceptible commuters to the working population at Shinjuku station in (A, B) and to the working population at Tokyo station in (C, D).
Figure S6: The 𝑹𝟎-centrality for every commuting pathway and residential station in the Tokyo metropolitan area, after vaccinating/quarantining 1,900 individuals from the largest working population at Shinjuku station. The 𝑅!-centrality for each commuting population and non-commuting population after vaccinating/quarantining 1,900 individuals from the largest working population at Shinjuku station, are given in accordance with the relation to the working population at Shinjuku station (currently the
(A) Commuting population (B) Non-commuting population 1 Home
Work Resident
1 Home
Work Resident
1 Home
Work Resident
1 Home
Work Resident
1 Home
Work Resident
(A-1) (A-2) (A-3) (B-1) (B-2)
Shinjuku
Commute to Shinjuku (2nd largest)
Commute to Tokyo (largest)
2 Home
Work Resident
2 Home
Work Resident
2 Home
Work Resident
2 Home
Work Resident
2 Home
Work Resident
(C) Commuting population (D) Non-commuting population
(C-1) (C-2) (C-3) (D-1) (D-2)
Tokyo
Commute to Tokyo (largest)
Commute to Shinjuku (2nd largest)
second largest susceptible work population after vaccination) and Tokyo station (currently the largest susceptible work population after vaccination) are given in (A for commuting population, B for non-commuting population) and (C for commuting population, D for non-commuting population), respectively. The schematic illustration above each panel describes its relationship. The color of dots indicates the number of susceptible commuters to the working population at Shinjuku station in (A, B) and to the working population at Tokyo station in (C, D).
Figure S7: The effect of countermeasures on the local final size of epidemic at each major station in the Tokyo metropolitan area. The change in the local final size of epidemic when the vaccination/quarantine is independently applied to the working population of each major station, Shinjuku, Tokyo, and Shibuya are given in (A), (B), and (C), respectively. Here, the vaccination/quarantine is applied to the relevant population only and the other populations are kept untouched. The result for vaccinating 0 (red circle dot), 1,000 (green circle dot) and 2,000 (blue circle dot) individuals are given and the local final sizes of epidemic at each work population are plotted against
(A1)
(A2)
(A3)
(B1)
(B2)
(B3)
(C1)
(C2)
(C3)
β =0.0001 β =0.0002 β =0.0008
Shinjuku Tokyo Shibuya
number of vaccinated/quarantined 0 number of vaccinated/quarantined 1000 number of vaccinated/quarantined 2000
its local population size. Each panel corresponds to results for different infection rate.
The sigmoidal profiles observed in Figure S1B are also evident here; for larger infection rate the transition point will shift to the smaller side. The overall shapes are not altered by the vaccination/quarantine, except for the relevant vaccinated/quarantined population.
This is because the number of vaccinated/quarantined is minimal compare to the total population size (i.e., less than 1%), so the effect of vaccination/quarantine is limited to the particular population only. A black arrow denotes the decrease of local final size of epidemic at each vaccinated/quarantined population.
