# <font color="darkblue">On Epidemic Indicators and Forecasting the New Normal: <br> Is measuring $R_0$ sufficient?</font>

by Christopher P. Monterola

---

In the appendix of Section 2, we have listed five reasons why $R_0$ is an important parameter in epidemic modeling. However, for purpose of developing actionable insights it is a very limiting measure and we will explain here why. We will also demonstrate how it can be modified in combination with other parameters as guide in lifting, imposing, or modifying quarantine protocols. 



## I. $R_0$ is not enough to ascertain outbreak, prevalence and number of susceptible population are equally critical.

The probability of extinction ($P_{ext}$), defined as the chance that an outbreak will persist in a SIR model of **homogeneous population** is derived in 1955 by [Whittle](https://www.jstor.org/stable/2333427?seq=1) as:

$$ P_{ext}(R_0, I_t) \sim   \left\{
\begin{array}{ll}
      (\frac{1}{R_0})^{I_t} & R_0 \geq 1 \\
        1 & R_0 < 1\\
\end{array} 
\right.  $$

where $I_t$ is the infected population at current time $t$. The chance of outbreak is just the opposite spectrum and is given by:

\begin{equation}
P_{outbreak}=1 - P_{ext} 
\end{equation}

or 

$$ P_{outbreak}(R_0, I_t) \sim   \left\{
\begin{array}{ll}
      1- (\frac{1}{R_0})^{I_t} & R_0 \geq 1 \\
        0 & R_0 < 1\\
\end{array} 
\right.  $$

**Four lessons from the above relation: **

1. $I_t$ is very critical. See simulation below.
2. $R_0$ is proportional to the effective transmissivity $\beta$, a contact parameter proportional to the number of susceptible population $N_s$ and/or density of $N_s$. The measured $R_0$ is derived with an assumed $N_s$.  



Note that if $R_0 >1$, even a few infections will force the system to an outbreak with exponential growth at the onset that's given by $I(t) = I(t_0)$ $\exp$$^{(\beta -\gamma)}$$^{(t - t_0)}$  where $R_0 = \beta /\gamma$ for the SIR model. See previous notebook for discussion. This leads us to the third point:

3. It is indeed critical to make sure that $R_0 < 1.0$ since this will guarantee that the highest number of infection will no longer exceed the numbers reported daily. However, consider always go back to lesson of item 2 since $R_0$ is a function of the current environment when it was measured and it can exceed 1.0 immediately after opening. 

4. If one wants to modify the quarantine when $R_0 > 1.0$, here is the most important policy measure that must be defined:   **Setting of an acceptable $P_{\texttt{outbreak}}$.**  The threshold number of cases ($I_{\texttt{threshold}}$ or the maximum cases reported that day) that will allow the extinction to go below the set $P_{\texttt{outbreak}}$ is given by:

\begin{equation}
I_{\texttt{threshold}} \sim \frac{- \log (1-P_{\texttt{outbreak}})}{log R_0}
\end{equation}

If $I_{\texttt{reported}}(t)$ at time $t$ is lower than $I_{\texttt{threshold}}$, then you satisfy your $P_{\texttt{outbreak}}$ criteria.




To better explain the above points consider the figure adapted from Allen et al (Letters in BioMathematics 2017, https://journals.plos.org/plospathogens/article?id=10.1371/journal.ppat.1003277). When $R_0 < 1$, $P_{outbreak} \sim 0$ while when $R_0 = 2$ with $I_t = 2$, $P_{outbreak} = 0.25$. 





### Policy Implication/Suggestion:

We need to work out simple rules that are implementable and explainable to our general constituents. Here is a suggestion if outbreak probability will be used as main indicator:

**Consider city/provincial level data, do not wait for probability of outbreak to be near zero in the national scale.** Make sure $R_t<1$ but after satisfying this condition note that the number of confirmed cases becomes a more important indicator than $R_t$ as $R_t$ is an unknown once you modify quarantine protocol. A prudent compromise on safety and accuracy is to assume $R_t \sim 1.2 - 1.4$ will be the immediate $R_t$ after opening, this is an acceptable estimate pre-quarantine for COVID19 with the assumption we better value social distancing after quarantine. 



```python
%matplotlib inline
from ipywidgets import interactive
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.collections import EventCollection
```


```python
def f(i_t, max_R0):
    
    plt.figure(figsize=(10,5))
    x = np.linspace(1, max_R0, num=100)
    plt.plot(x, 1-np.power(1/x,i_t), 'o-')
    plt.ylim(0, 1.1)
    plt.xlabel("R0",  size=14)
    plt.ylabel("Probability of Outbreak", size=14)
    plt.show()

interactive_plot = interactive(f, i_t=(0, 50,1), max_R0=(1.1,5,.1))
output = interactive_plot.children[-1]
output.layout.height = '350px'
interactive_plot
```


    A Jupyter Widget



```python
def f(P_outbreak, max_R0):
    
    plt.figure(figsize=(10,5))
    x = np.linspace(1.01, max_R0, num=100)
    plt.plot(x, -np.log(1-P_outbreak)/np.log(x), 'o-')
    plt.xlabel("R0",  size=14)
    plt.ylabel("$I_{threshold}$", size=14)
    plt.show()

interactive_plot = interactive(f, P_outbreak=(0.05, 1, .05), max_R0=(1.0,3,.1))
output = interactive_plot.children[-1]
output.layout.height = '350px'
interactive_plot
```


    A Jupyter Widget


### Illustration 1.

If the government sets the acceptable $P_{outbreak}$ to $49\%$ in modifying the quarantine, if $R_0 = 1.1$ you need to have cases to be lower than $I_{threshold} \sim \frac{- \log (1-P_{outbreak})}{log R_0} = 7$ cases. 

Important caveat: The model approximates a homogeneous population. See Section 4 below for heterogenerous population which increases the $I_{threshold}$, i.e. make it less stricter. 


```python
print(-np.log(1-0.49)/np.log(1.1))
```

    7.064770569588317


---

## II. Do not be obsessed in interpreting $R_0$.  

We are often obsessed with the idea of $R_0$, whether or not it is less than one. We should not be obsessed by these fluctions. Let's illustrate this point by answering this question: 

> If we observe $R_0$ to be below 0.5, can we now assume it is safe to lift the quarantine as compared to 0.76?

No. The most simple argument is to use actual data from the Philippines' National Capital Region (BCR). Let's use April 08 as reference point when we had the lowest national $R_t$ to date. [Interact with the computations here](http://www.complexity.asia/covid-rt/) (designed for viewing in desktop/laptop).


Notice that four days later, the $R_t =0.11$ in San Juan eventually increased to $1.56$ on April 12, while $R_t =0.86$ in Las Piñas lowered to 1.47 on April 12. See table below. 


| City |  April 08 ($R_t$) | April 12 ($R_t$) |
| --- | --- | --- |
|San Juan | 0.11 | 1.56 |
|Pasay | 0.42 | 1.55 |
|Makati | 0.49 | 1.00 |
|Manila| 0.63 | 0.89 |
|Las Pinas | 0.86 | 1.47 |


Is this a general result? No. However, note that cases with $R_t <1$ are non-trivial; and, it is seemingly futile to segmentize since the dynamics will significantly change once quarantine is lifted. 

For practical reasons, it cannot tell us much except that it is way better than when $R_t > 1$ -- the theoretical threshold for an epidemic model. The transition from zero ($R_t <1$) to $P_{ext} \sim (\frac{1}{R_0})^{I_t}$ ($R_t <1$) is the most scientifically valid policy point. 

Another thing to note is the high correlation of the reported $R_t$ among the NCR cities. In general, the trend of NCR is leading the trend of cases outside NCR. This is an expected result since NCR is the local "epicenter" of COVID-19. Consequently, NCR cities are less correlated than others. This observation can be explained as a combination of: 

(1) system riding on the same testing trend,    
(2) proximity of cities and hence higher chance of population mixing,  and   
(3) similarity of quarantine dynamics imposed. 

You will find a lot of scientific articles on interpreting $R_t <1$ but most of those articles are based on a specific data and not generalizable. Our empirical data speaks to us and this is the most prudent way of interpreting $R_t <1$.



---

## III. Homogeneous population assumption is a baseline approximation and in terms of spreading likelihood assumes the worst case of real-world heterogeneous population 

In general, the population is not homogeneous, **how do we estimate then the impact of heterogeneity?** 

We can approximate it using the curve-fitting result of [Hartfield *et al*](https://journals.plos.org/plospathogens/article?id=10.1371/journal.ppat.1003277) using the idea of [Smith *et al*](https://www.nature.com/articles/nature04153) as a function of different heterogeneity factor $k$. The homogeneous case is when $k =1 $; and, as $k$ decreases to zero, the system becomes more heterogeneous.  As a result of heterogeneity of population (spatial location, temporal location,  daily routine), the contact dynamics of people becomes more constrained. What do we mean by this? If the system is homogeneous, we can expect that the likelihood of one person meeting or encountering another  person at random is the same as meeting any other random person or a friend, a relative, or someone in the family. However, we know that this is not the case. 

Hartfield *et al* estimated that  $I_{threshold} \sim \frac{- \log (1-P_{outbreak})}{log R_0}$ should then be updated as:

\begin{equation}
I_{threshold} \sim \frac{- \log (1-P_{outbreak})}{\log (R_0)} (.334 + \frac{.689}{k} + \frac{.408}{R_0} - \frac{.507}{k R_0} - \frac{.356}{R_0^2} + \frac{.467}{k R_0^2})
\end{equation}

The result is shown in the figures below. Notice that heterogeneity relaxes a bit the condition on the number of confirmed cases that need to be achieved. 


### Illustration 2.

Refer to Illustration 1 above and use $k=0.16$ (heterogeneity of SARS outbreak in Singapore 2003). Again, if the government sets the acceptable $P_{outbreak}$ to $49\%$ in modifying the quarantine with $R_0 = 1.1$, instead of lower than 7 cases for homogeneous, our new threshold for heterogeneous population is 30 cases (see computation below). 



```python
P_outbreak =0.49
k=.16
x = 1.1 #R0

I_threshold = -np.log(1-P_outbreak)/np.log(x)
I_threshold =I_threshold * (.334 + (.689/k) + (.408/x) - (.507/(k*x)) - (0.356/x**2) + (.467/(k*x**2)))
print(I_threshold)
```

    30.014342042327574


### Policy Implication/Suggestion

If you want to prioritize safety of population due to COVID19 infections, use the homogeneous case as your baseline. If you want to include economic constraints and relaxed a bit your safety anchoring use k ~ 0.2 to 0.4, a heterogeneity factor that covers the range of SARS to Spanish Flu to Smallpox.

Find in the next two simulations how $k$ changes the outbreak probability.


```python
def f(P_outbreak,k):
    plt.figure(figsize=(10,5))
    x = np.linspace(1.01, 3, num=100)
    I_threshold = -np.log(1-P_outbreak)/np.log(x)
    I_threshold =I_threshold * (.334 + (.689/k) + (.408/x) - (.507/(k*x)) - (0.356/x**2) + (.467/(k*x**2)))
    
    if k==.16:
        print('SARS, Singapore 2003')
    if k==.25:
        print('Measles, US 1997-1999')
    if k==.35:
        print('Smallpox, Europe 1958-1973')
    if k==.65:
        print('Smallpox, England 1966')
    if k==1.0:
        print('Homogeneous population')
    plt.semilogy(x, I_threshold , 'o-')

    plt.ylim(1, 1000)
    plt.xlabel("R0",  size=14)
    plt.ylabel("$I_{threshold}$", size=14)
    plt.show()

interactive_plot = interactive(f, P_outbreak=(0.05, 1, .05), k =(0.05, 1, .01) )
output = interactive_plot.children[-1]
output.layout.height = '350px'
interactive_plot
```


    A Jupyter Widget



```python
def threshold(k, P_outbreak, x):
    I_threshold = -np.log(1-P_outbreak)/np.log(x)
    I_threshold =I_threshold * (.334 + (.689/k) + (.408/x) - (.507/(k*x)) - (0.356/x**2) + (.467/(k*x**2)))
    return I_threshold

def f(P_outbreak,max_R0):
    plt.figure(figsize=(10,5))
    x = np.linspace(1.01, max_R0, num=100)
    
    plt.semilogy(x, threshold(0.16, P_outbreak,x), 'r', label='k =.16, SARS, Singapore 2003')
    plt.semilogy(x, threshold(0.25, P_outbreak,x), 'g', label='k =.25, Measles, US 1997-1999')
    plt.semilogy(x, threshold(0.35, P_outbreak,x), 'b', label='k =.35, Smallpox, Europe 1958-1973')
    plt.semilogy(x, threshold(0.65, P_outbreak,x), 'm', label='k =.65, Smallpox, England 1966')
    plt.semilogy(x, threshold(1.0, P_outbreak,x), '-', label='k =1.0, Homogenous case')
   # plt.semilogy(x, threshold(5.0, P_outbreak,x), '--', label='k =5.0, Homogenous cramped case')

    plt.legend()
    plt.ylim(.5, 1000)
    plt.xlabel("R0",  size=14)
    plt.ylabel("$I_{threshold}$", size=14)
    plt.show()

interactive_plot = interactive(f, P_outbreak=(0.05, 1, .05), max_R0=(1.5,5,.1))
output = interactive_plot.children[-1]
output.layout.height = '350px'
interactive_plot
```


    A Jupyter Widget


## IV. Understanding the contact network driven by urban mobility after quarantine is lifted is the key in strategizing in this new normal. 

Together with mass testing, we can further improve our strategy.

Some of our ideas are previously published here: 

(1) [Importance of small-degree nodes in assortative networks with degree-weight correlations](http://www.complexity.asia/covid-rt/docs/network-epi/papers/pre-2017/)
Ma, S., Feng, L., Monterola, C. P., & Lai, C. H. (2017). Importance of small-degree nodes in assortative networks with degree-weight correlations. Physical Review E, 96(4). doi:10.1103/physreve.96.042308




(2) [Encounter network in disease spreading](http://www.complexity.asia/covid-rt/docs/network-epi/papers/proc-cs-2015/)
Ramli, M. A., & Monterola, C. P. (2015). The Resilience of the Encounter Network of Commuters for a Metropolitan Public Bus System. Procedia Computer Science, 51, pages 2117-2126. doi: 10.1016/j.procs.2015.05.482



However, the use of the above science in the Philippine setting entails: 

1. Mapping the contact network in our transport system,    
2. Mapping the contact network for various ages or age groups,    
3. Or, in general, mapping the contact network of various segments of society.   
    
Here are some approximations on how we should estimate the limit for mass transport:


A. **Transportation limit.** Maintaining social distance in previously cramped MRT means reducing the capacity by $1/3.75$ for a 1 meter social distancing protocal (from initially ~ 3 by 5 = 15 you reduce the space to 4/15 =1/3.75) and 1/12.5 for a 2 meter social distancing protocal (5 by 10 = 50 you shrink the space by 4/50 =1/12.5). This does not take into account the queue of the MRT in usual days. A good compromise is probably to initially set the transport limit to 1/15 to 1/20 of the commuters. This means that we should have a continuing workforce calculated under modified ECQ in the range of 5 to 6.7%.

B. **Office space limit.**  We expect that transport sets the lower limit (MRT is more cramped than office), so this can be relaxed. However, regular testing should be be done to the "super-spreaders" who can be teachers in classrooms, cashiers in groceries, or errand staff from offices. Our recommendation, every 3-4 days to capture the early infectious state.




### Policy Implications/Suggestions

1. If quarantine is lifted, start the modification by limiting the workforce in NCR by 5 to 7% of workers per office. This can be more relaxed in other cities or provinces having a more efficient transport system per capita. Priotizing employees residing near the office will have the benefit of relaxing further transport congestion and spreading of virus due to prolonged contacts.

2. if quarantine is lifted, test potential super-spreaders every 3-4 days, (maybe) this cost can be shouldered by the institution who would like to fully return to operation.





```python

```

6084e9d7102f720c763f490fde2c841e45697248
ipynb
Jupyter Notebook
notebooks/Part IV Epidemic Indicators - Insights and Practical Considerations .ipynb
complexityasia/epimod
9f14113916f786fd97bd40a0eac0288584b4d12b
MIT
notebooks/Part IV Epidemic Indicators - Insights and Practical Considerations .ipynb
complexityasia/epimod
9f14113916f786fd97bd40a0eac0288584b4d12b
MIT
notebooks/Part IV Epidemic Indicators - Insights and Practical Considerations .ipynb
complexityasia/epimod
9f14113916f786fd97bd40a0eac0288584b4d12b
MIT
Qwen/Qwen-72B
1. YES
2. YES