Effects of control measures and their importance on COVID-19 transmission dynamics

For more than a year, governments around the world have attempted to control the COVID-19 pandemic. Control measures such as social distancing, face mask wearing, business/school closure, city or transportation lockdown, ban of mass gathering, population education and engagement, contact tracing, and improved mass testing protocols are being used to contain the pandemic. Currently, there are no studies to date that rank the importance of these measures so that the governments may allocate and target their resources towards the most effective control measures. In this paper, we propose a Discrete Time Markov Chain model that captures the above control measures and ranks them. We also show that the importance of the measures change overtime and depends on the stage of the transmission dynamics, as well as the environment. For example, contract tracing is known to be a powerful measure to effectively control the pandemic, however its influence is dynamic in nature. Our results show that contact tracing is indeed helpful during the early stage of the pandemic, but becomes less important after a vaccination program takes effect. If implemented, our novel and unique model may assist many countries in their crucial pandemic control decisions.


Introduction
Different countries use different control strategies to combat the COVID-19 pandemic. China and New Zealand implemented a seal and eradicate strategy, where they applied strict measures at the very early stage to stop local transmission fast, and then sealed their borders to protect the epidemic from reoccurring. The US and UK choose to balance between controlling the pandemic and its economic impact by implementing various measures to keep the effective reproduction rate Rt low, while attempting to end the pandemic with their advanced vaccination programs. Widely accepted control measures used in pandemic control are social distancing, face mask wearing, business/school closure, city or transportation lockdown, ban of mass gathering, population education and engagement, contract tracing, vaccination programs, viral testing advancements, and targeted mass testing. In order to formulate an effective control strategy, decision makers need to understand the impact of each control measure to the transmission dynamics and be able to identify the control measures that should put more focus on and allocate more resources into.
Huang et al. 1 presents a SEIR compartmental model to show that the quarantine and isolation measures, as well as testing delays, could cause the largest difference in daily confirmed cases. Zhang et al. 2 found that limiting the mobility of the general population is more important than early detection and isolation during early stages of the outbreak in mainland China. Liu et al. 3 found that earlier and stricter lockdown implementation yields better control than later or more relaxed measures. Most studies focus on specific control measures and show their positive benefit. So far, there has not been a study that has attempted to aggregate all of these measures into a single model to see how they are interrelated and to compare the levels of their benefits. In this paper, we propose a Discrete Time Markov Chain model that captures most-if not all-widely accepted control measures and rank them so the policy maker can appropriately allocate their resources to the measures that are most beneficial.

Model
Some studies 4,5 use Markov chain with four states, namely Healthy, Active, Recovered and Death, to model COVID-19 evolution. Although having few states makes the model easier to manage, these models fail to accurately capture important dynamics of the epidemic, such as the chance of viral transmission and chance of recovering, which are nonstationary in nature and depend on how many days the patient has been infected. Our proposed model is more inclusive by having a greater number of states to capture the temporal information on how many days since a patient became infected, allowing nonstationary representation of transmission probability and recovery probability.
Our model can also capture the nature of population where there may be some portion of the population that will not cooperate with the government-mandated control measures to prevent viral transmission. Hence, population is classified into two classes, which are "cooperative" and "uncooperative". Let a control measure ko be the proportion of cooperative people in the population. With effective communication and promotional campaigns, the government might be able to improve ko.
There are a significant portion of asymptomatic cases. Most asymptomatic cases will never be detected but can still spread the virus similar to symptomatic patients. Moreover, those with mild symptoms that have recovered by themselves may not get tested, also contributing to the undetected rate in the population. These two types are classified as "undetected infected" class. Additionally, there may be some detected and infected person, possibly with mild symptoms, that did not comply with self-isolation. Those who did not self-isolated are also assigned to the undetected infected class as well since they may possibly go out and spread the virus in a similar way as the undetected infected individuals. Thus, the undetected infected class include asymptomatic and mild symptom cases who do not get tested, and those that tested positive but failed to isolate themselves. Recent studies have estimated asymptomatic cases to be around 20% [6] to at least 33% [7] with wide range from 4% to 52% [8], 43-43.2% [9] and [10], and 21.9%-35.8% [11]. Without loss of generality, we set the asymptomatic proportion to be 35% and implement a control parameter kd as a percentage adjustment to the asymptomatic proportion, in order to represent an ability to increase asymptomatic and symptomatic cases detection and isolation. That is, the proportion of undetected infected over total infected is 35% -kd. For example, if we can detect 30% of asymptomatic cases and can isolate 80% of symptomatic case, we have kd = 30% (35%) -20% (65%) = -2.5%. We estimate that kd can range from -10% to 30%, depending on the effectiveness of the control measures being implemented.
Our Discrete Time Markov Chain model at time t has the following states. SCt = Susceptible and cooperative with control policy on day t SUt = Susceptible and uncooperative with control policy on day t UC i t = Undetected infected for i days and cooperative with control policy on day t UU i t = Undetected infected for i days and uncooperative with control policy on day t TC i t = Detected or will be detected infected for i days and cooperative with control policy on day t TU i t = Detected or will be detected infected for i days and uncooperative with control policy on day t Mt = Cumulative immune on day t Vt = Cumulative vaccinated on day t Dt = Cumulative dead on day t, where i = 1, …, I , I is the maximum number of days 1 that an infected person can transmit the virus. In other words, fu(i) and fd(i) should be negligible when i > I.
The probability for an infected individual to spread the virus after he/she became infected for i days follows temporal transmission probability mass functions, fu(i) and fd(i), for undetected and detected infected individuals, respectively. Let R0 be the basic reproduction number, which is the expected number of cases directly generated by one infected individual. Let us define an aggregated epidemic control measure kr(t) as the measures implemented on day t that directly affects the reproduction number, such as mask wearing, social distancing, and lockdown. For example, Oraby 5 found that lockdown reduces R0 by 64-85% among 155 countries. By design, only the cooperative population will implement this kr(t) measure while the uncooperative population will not. Hence, the expected number of new cases generated from the contagious population on day t can be computed as Et = (Si UC i t fu(i) + Si TC i t fd(i)) · R0 kr(t) + (Si UU i t fu(i) + Si TU i t fd(i)) · R0.
Undetected cases are asymptomatic or mild symptomatic therefore we assume that they will all recover and become immune. For detected cases, some of them will die. Recent global case fatality rate (CFR) 12 is 2.71%. Hence, the probability that patients in the detected group will die is assumed to be 2.71%, while the remaining 97.29% of them will recover. The estimated CFR of 2.71% is also in a reasonable range with other studies 2 . As a result, the 1 In the Methods section we show that I = 21 2 Wu et al. 24 found that for those that have symptoms, 81% have mild to moderate symptoms, while 14% have severe symptoms and 5% are critically ill. Assuming 10%-25% of severe symptoms and 100% of critical symptoms are in ICU, this 2.71% CFR translates in to immune state Mt+1 equals Mt +(UC I t + UU I t)+ 97.29% (TC I t + TU I t) and the death state Dt+1 equals Dt + 2.71% (TC I t + TU I t).
One important control measure is contact tracing. We define control parameters kl as tracing delay, the time between isolation of an index case and isolation of its contacts, and kc as contact tracing coverage, the proportion of contacts detected and isolated, which is equal to the percentage of cases that conduct contact tracing multiplied by the proportion of contacts that can be isolated over all traced contacts. Lastly, one of the most effective control measures is vaccination. Let kv(t) be the vaccination coverage being administered up to time t, which is equal to the percentage of population who received vaccination up to time t multiplied with vaccine efficacy, e.g. 95% for Moderna 13 and Pfizer-BioNTech 14 and 72% for Johnson & Johnson 15 . More details of contact tracing and vaccination modeling can be found in Methods section.

Case study I: USA as of March 21, 2021
As an example, we apply our model with recent data from the United States. As of March 21, 2021, there are 81,415,769 people 16 who received a COVID-19 vaccination. Since the current US population 17 is 331,002,651, the percentage of the population that has received at least one dose of the vaccine is 24.6%. At this time, the United States vaccination capacity follows an exponential growth rate with last 10 days average daily growth rate of 0.5%, derived from the data from [18]. According to a survey 18 , 47.2% of people are willing to take a COVID-19 vaccination if it was available to them. Hence, we should not expect kv(t) to go near 1.0 easily. Therefore, we assume for a base case that the vaccination program will continue at 0.5% growth rate until it reaches 60% of population. For the worst case, the program will continue at 0.5% growth rate until it reaches 45% of population. For the best case, the program will continue at 0.55% growth rate until it reaches 75% of population. Since more than 98% of vaccines administered 18 in the US are Moderna or Pfizer/BioNTech with efficacy 95% so we will assume vaccine efficacy of 95%.
Gramlich's 2018 survey 19 found that 75% of Americans would cooperate with each other in a crisis. This 75% number is also consistent with survey 20 of mask use in March 2021. Hence, we use 75% as a base case for ko and 50% and 95% for worst and best cases, respectively.
From the current CDC recommendations 21 for the purpose of advance public health planning, the best estimate of R0 to use in our model is 2.5. We then fit the parameters with actual data in order to obtain base case of control measure kr(t). By minimizing sum of square error in fitting the 7-day moving average of new case numbers from actual and model for the 42.34% -31.88% ICU death rate (=2.71/(5+0.1*14), which is in line with ICU death rate from previous studies. For less crowded health systems, Chen et al. 25 shows 22% death rate among ICU patients and Gold et al. 26 provides ICU death rate in a range of 37%-48.7%. In a crowded health system, however, the ICU death rate can be as high as 78% [27]. past 3 weeks, we found that stationary kr(t) at 0.4354 gives the best fit. As a result, our base case for kr(t) is 0.4354. Table 1 show a summary of control parameters and their assumption.

Table 1: Summary of control measure parameters and their assumptions
The United States is quite advanced in their vaccination program because as of March 21, 2021, nearly 25% of its population have received at least one dose of the vaccine. Figure  1 shows the number of daily new cases, under the base, worst, and best cases of control measure kr. In the base case when kr = 0.4354, it will take 136 days to get zero new cases. If a lockdown is imposed where kr becomes 0.2177 then it will shorten this duration to become 95 days. When kr is at 0.8708, e.g. people meeting at doubling current rate, the number of new cases transitions to an upward trend but later descends in mid-May when the percentage of the population vaccinated increases, yielding a total of 275 days to get zero new cases. Extend Data Fig. 1 shows the number of daily new cases, under different vaccination program scenarios. All three scenarios would not yield any significant difference in terms of the daily new case numbers. At 0.55% daily growth rate and max at 75%, it would take 124 days to have zero new cases, while at 0.5% daily growth rate and max at 45%, the zero new case duration will increase to 160 days. If the vaccination program stops on March 22, 2021, the daily new cases will still go down if the control measure level is maintained at the same level as in March 2021. This is because the effective reproduction rate Re is already less than 1.0 from the 25% vaccination coverage. The zero new case duration, however, will extend to 289 days. There is a herd immunity threshold 22 which is required for control of transmission. However, as seen from the result, the benefit of vaccination can be attained right away without having to wait to reach any specific threshold, as long as the control measures to keep Rt down are still enforced.  Fig. 2 is a tornado diagram displaying the number of days to zero new cases from the best and worst cases of each control measure. For example, reducing kr(t) to 0.2177 (best case) results in 95 days to zero new cases while increasing kr(t) to 0.8708 (worst case) results in 275 days to zero new cases. The diagram shows that, in the US where its vaccination program is quite advanced, to further improve the speed of pandemic control, the policy makers should focus on 1) lower the effective reproduction even more or at least keep it at the base level by effectively inform the people to maintain social distancing and wearing face mask habit, 2) improve the cooperation from the population by education and engagement, and 3) reducing testing delay. On the contrary, the contact tracing program and asymptomatic detection ability now become less important and the control speed would not be affected much, e.g. just only 11 days, by improving their parameters.

Case study II: Thailand as of January 6, 2021
To understand the effect of the control measures in a different environment, we apply the model with parameters estimated from Thailand data as of January 6, 2021. Table 2 represents worst, base, and best case scenarios for the analysis. From fitting with the actual data and set kr(t) at 1.0, we found that R0 for the 2 nd wave in Thailand was 2.026. Thus, the worst case for kr(t) is 1.0 which is the same level during the zero case period before the start of the second wave. Rotejanaprasert et al. 23 estimated Rt in Bangkok and shows that it can go down as low as 0.5 in May 1, 2020. This implies kr(t) of 0.2, which is same level when Thailand was in lockdown during April/May 2020. For the base case, we assume kr(t) in the middle, which is 0.6. Fig. 3 shows daily new case numbers when similar lockdown measures implemented in April 2020 (kr(t) = 0.2) were re-implemented after the 2 nd wave was detected for 2, 9, and 16 days delay. The graph shows that the longer the control measure was delayed, the harder and longer it takes to bring the new cases down. If the lockdown was implemented in 2, 9, 16, and 25 days, it would take 30, 43, 55, and 67 days to bring the new cases to zero, respectively.   As seen from Fig. 4, even though Thailand vaccination program lags behind many other nations, if it implements a medium control level kr at 0.6, it will begin to see a positive effect of the vaccination program when kv reaches 5%, corresponding to 4.6 million people received vaccination. At this point, the number of daily new cases will start to descend. This is because the effective reproduction rate is near 1 so just only a small reduction in susceptible population would be good enough to bring the effective reproduction rate to be below 1.

Discussion
The tornado diagram of Thailand in Fig. 5 indicates that, even though Thailand currently has low number of new cases, it will take more than a year to get zero new cases if all control measures are still kept at their base levels. Some serious efforts are needed to be put on to improve some -if not all -of these control measures to better levels. Most benefits would be obtained from implementing measures that reduce reproduction rate and/or testing delay, such as 1) ban of unessential mass gathering, 2) increasing mobile testing units that offer free COVID-19 test to anyone with symptoms or who live/work in a highly populated area, 3) quickly, practically, and effectively isolate those with positive test results so that they can no longer spread the virus. In addition, contact tracing and asymptomatic detection ability are still needed to be implemented and improved. Moving any of their parameters from its worst case to its best case will affect the number of days to zero cases by more than 100 days. Most importantly, there are three measures; which are reducing reproduction rate, increasing cooperation level, and vaccination program; that need special attention to make sure that they do not fall to their worst-case levels otherwise the pandemic can be out of control.

Temporal transmission probability of undetected cases
The undetected group includes asymptomatic cases, cases with mild symptoms but untested, and cases with mild symptoms and positive test results but failed to isolate themselves. Zou 28 found that the viral load detected in the asymptomatic patient was similar to that in the symptomatic patients. We assume they have the same generation time profile as symptomatic cases and use Weibull distribution estimated from Ferretti et al. 29 . They used maximum likelihood estimations to infer the generation time distribution and fit it with a Weibull distribution with a mean of 5.5 days and standard deviation of 1.8 days. This translates to a Weibull distribution with a shape parameter of 1.6749 and scale parameter of 6.1577. We apply this distribution to R0 to assign the probability that one infected person can spread the virus to other persons in each day, as shown in

Temporal transmission probability of detected cases
For detected cases, to estimate the probability that an infected person can spread the virus in each day, we combine incubation period distribution 3 estimated from Lauer et al. 30 with the infectious profile 4 estimated from He et al. 31 and updated by Ashcroft et al. 32 in order to deduce generation time, as shown in Table 1. The rows in the table represent possible cases with different incubation periods, that can range from 1 to 21 days, with associated probabilities from Lauer et al. 30 shown in the first column. Each row illustrates the probability that an infected person with the given incubation period can transmit the virus on each day during the course of the infection. The infectiousness profile data are from He et al. 31 and updated by Ashcroft et al. 32 . The "weighted average" row shows the incubation probability weighted average of the infectiousness probabilities. The normalized numbers are shown in the "normalized" row in order to make the probabilities sum to one. Note that we cut of the tail of infectiousness profile after 21 days because, when multiplying with incubation period probability, the results are negligible. Hence, we assume that the maximum number of days that an infected person can transmit the virus, I, is 21. If an infected person does not recover in 21 days, he/she should be hospitalized and should no longer be able to transmit the virus. In many cases, we can improve the control of pandemic by imposing isolation after detection. This is usually done by quarantine, hospitalization, or self-isolation. Let us define the time lag between the time an infected person develops symptoms and the time that he/she is isolated as testing delay or kg, which can be vary from 0 to 7 days 34 , while the infectiousness probabilities after kg day from the onset will be set to zero. The example truncated probabilities with kg = 3 are shown in Table 2. The truncated probability numbers then are normalized with the factor from no isolation cases since isolation after detection should have no impact to the transmission probability before isolation.
after symptom onset would be safe to not be in isolation. Some parameters of the infectious profile were updated in Ashcroft   The results of the cases with 1 -3 days isolation after symptom onset are shown in Fig.1

Contract tracing modeling
Let us define incubation period distribution, which is estimated from Lauer et al. 30 , as P, where p(j) represents the probability that the incubation period is j days. Suppose a detected infected individual has an incubation period of j days, the time that he/she will be isolated is j+kg days, and the time that proportion kc of his/her contacts will be isolated is j+kg+kl days. On a given day t, the expected numbers of patients with j days incubation period that will be detected and isolated on that day are ( ) # .$/ 0 for cooperative group  7 8 − 8 + < = + = + ( ) # .$/ 0 2 4 ( ). It will take kl days until their contacts are isolated. As a result, on day t, contact tracing will reduce the number of contagious people in each state of the Markov chain as follow.
Suppose that a contact of an index case was isolated on day t, this means that the index case was detected and isolated on day t-kl. On day t-kl, the index case has got the virus for j+kg days, where j is his/her incubation period. The first day that he/she got the virus would be day t -kl -( j+kg -1). Since then he/she has been potentially spreading the virus for i days, i = 1, …., j+kg-1, corresponding to day t -kl -( j+kg -i), from day t -kl -( j+kg -1) to day t -kl -( j+kg -(j+kg-1)) = t -kl -1. Hence, on day t, the contact who got the virus from the index case on the index case's i th day of infection, which was day t -kl -( j+kg -i), would be infected for t -(t -kl -( j+kg -i)) = kl + j + kg -i days. This explains the superscript of UC in equation (1).
We consider only the cases when the contact has been infected less than I days, i.e., kl + j + kg -i ≤ I, because after I days the contact can no longer transmit the virus. Thus, i ≥ kl + j + kg -I. We have that, i ≥ max(1, kl + j + kg -I ). Moreover, the index cases can spread the virus at most I days. Hence, i ≤ min(I, j+kg -1). This explains the range of i in (1).
The reduction in the numbers in states UUt follows the same form as UCt in (1), except that SC in (1) is replaced by SU. The reductions in numbers in state TCt and TUt follows the same forms as those of state UCt and UUt, respectively, except the term (35%kd) is replaced by term (65% + kd) instead to appropriately represent detected group. Moreover, the sum of the reductions from state UCt and UUt are added to state Mt because these isolated contagious contacts will not be able to transmit the virus anymore. However, only 97.29% of the sum of the reductions from state TCt and TUt are added to state Mt since 2.71% of them will die, which in this case are added to state Dt.

Vaccination modeling
Due to possibility of reinfection, countries such as the UK 35 and US 36 and currently provide vaccination to their populations regardless of whether they had already been infected or not. We define µ to be the proportion of COVID-19 recovered people that can get vaccination. Based on the current vaccination policy, we assume that µ = 1. Let ke be the number of days for the vaccine to give full protection. CDC 37 currently suggests that ke be 14 days. As a result, the number of people receiving vaccination at time t is kv(t)(N-Dt) -kv(t-1) (N-Dt-1). The population that can receive vaccination at time t is N-Dt-Vt-(1-µ)Mt. Therefore, the number of cooperative susceptible persons SCt to get vaccine at time t is (kv(t)(N-Dt) -kv(t-1)(N-Dt-1))· SCt / (N-Dt-Vt-(1-µ)Mt). These vaccinated individuals would still be susceptible for the next ke days as of them may become infected before the vaccine reaches its maximum protection. We assume that, during the ke days before the vaccine takes effect, the chance that susceptible individuals to become infected is the same regardless of whether they took the vaccine or not. Hence, at time t, the number of cooperative susceptible SCt will be reduced by (kv(t-ke)(N-Dt-ke) -kv(t-ke-1)(N-Dt-ke-1))· SCt / (N-Dt-ke-Vt-ke-(1-µ)Mt-ke) and moved to Vt. The number of infected individual that took vaccination at time t-ke is (kv(tke)(N-Dt-ke) -kv(t-ke-1)(N-Dt-ke-1))· (SCt-ke -SCt) / (N-Dt-ke-Vt-ke-(1-µ)Mt-ke), which is subtracted from Mt and moved to Vt. The formula for state SUt follows the same fashion as SCt but with SU replacing all SC.