Three-level algorithm to control COVID-19 spreading

As the main ways of the coronavirus (COVID-19) spreading are air and physical contact, actions to mitigate and suppress its spreading must be developed in order to change population behavior, providing eﬃcient control strategies. Here, these actions are described as a simple heuristic framework to establish public policies. The possible eﬀects of these actions are modeled as a stratiﬁed control algorithm applied to an extension of the Susceptible-Infected-Removed (SIR) compartmental model. The model dynamics is analyzed and validated with simulations considering parameters obtained from epidemiological data from Brazil and Uruguay. Model simulations allow the proposition of a three-level control strategy that, as the numerical experiments show, presents promising results for the development of strategies to reduce the spread of COVID-19.


Introduction
In the beginning of 2020, the epidemic of a new disease changed the history of humanity, being responsible for a global health and economic crisis. COVID-19 is a new coronavirus progeny, namely SARS-CoV-2 and was firstly reported in Wuhan, China, on December 31, 2019 [1,2]. A few months later, a report published by the World Health Organization qualified the disease as a pandemic due to its worldwide prevalence [3].
Transmission of SARS-CoV-2 occurs predominantly through physical transportation of contaminated droplets of secretions from an infected individual to an uninfected person, although the role of aerosol transmission and contaminated surface contact is currently unknown [4]. It is a disease with lower lethality when compared to its high transmission capacity, with the spreading aggravated by the average incubation time and due to both symptomatic and asymptomatic people being able to transmit the disease [5].
The disease has several forms of manifestation and most patients have rapid resolution, many of whom are asymptomatic. The initial symptoms are fever, cough and breathing difficulties. However, some patients with comorbidities may develop complications [6], requiring hospitalization, intensive treatment and mechanical respiratory ventilators, provoking a severe crisis in the health care system caused by the excessive number of hospitalizations.
To avoid overcrowded hospitals, different interventions aimed to slow the rapid evolution of the pandemic have been adopted by many countries. Strategies such as isolation of infected patients, mask wearing, regular hand hygiene, social distance, quarantine and total blocking of areas [7,?,?,?] were adopted all over the world.
Mathematical models are used to simulate the transmission of coronavirus and are helpful to understand the dynamical behavior of an infection, providing tools to estimate the duration and peaks of infection outbreaks and to design efficient control strategies [11,12,13]. Epidemiological mathematical models to be used here are present in health literature [14,15] related to several diseases, being helpful to define public actions to reduce COVID [16,17,18,19,20,21,22].
Development of macroscopic compartmental models for disease spreading has been research field since the Kermack and McKendrick proposition [23,24,25], dividing population into groups such as in Susceptible -Infected -Removed (SIR) or in Susceptible -Exposed -Infected -Removed (SEIR) models. Several modifications of the SIR model [26,27,28] have been proposed for epidemic modeling [29,30,31,32], being useful to public health policy makers [33].
Due to the urgency of COVID-19, researchers around the world accepted the challenge of outlining mathematical models for the new epidemic. Several mathematical models have already been formulated for the population dynamics of COVID-19 in several countries [10,34,35,36,37,38,39,40] and pioneer methods are structured upon Machine Learning and statistical models such as decision trees and linear regression [41] or even more powerful ones like artificial neural networks [42,?], showing promising results. Machine learning has also been used as an assistance to early detection [44] and prediction of severity [45] of COVID-19.
However, this paper considers a more classical and established approach, using the adjustment of social contact rate characterized by the control of the β parameter to mitigate COVID-19 spread. For this, a three-level controller is implemented and its relationship to social distance is validated. An extension of the compartmental SIR model is considered, the Susceptible -Infected -Removed -Dead (SIRD) and the method is based on control theory, helping to define public policies to decrease the COVID-19 propagation [46].
Remaining paper is ordered as follows: in section 2, the equations and hypothesis of the model are presented. In section 3, the analysis of model dynamics is exhibited, the equilibrium points are determined, basic reproduction rate R 0 for the SIRD model is defined and analyzed with two simulations. Then, in section 4 the model is studied for the sake of establishing a baseline for the controller, with parameters derived using epidemic data from Brazil and Uruguay. Finally, section 5 outlines the controller and its effectiveness is investigated for a hypothetical scenario and the conclusions are shown at the end, in section 6. All simulations are performed with MATLAB Simulink.

Model formulation
The proposed model is a modification of the original SIR compartmental model [47], including a dead population compartment, here called Susceptible -Infected -Removed -Dead (SIRD) model, shown in figure (1). The model does not account for natural births and deaths as the Dead compartment being related only to deaths due to complications caused by the coronavirus infection. Model parameters are described in table (I). β is the average number of contacts that result in contamination, per unit of time; γ is the recovery rate, i.e. the rate at which individuals leave the infected compartment after being recovered. Consequently, γ −1 represents the recovering time [48]. Additionally, Ω is the fraction of the infected population that dies, per unit of time (mortality rate).

Parameter
Parameter description β Average contact rate.
Mean infectious period. Ω Mortality rate. As it can be observed, the model presents some simplifications: a recovered individual is supposed to acquire immunity and does not return to the susceptible compartment becoming "Removed" ("R" from SIRD); -parameters β, γ and Ω are considered to be constants, despite depending on individual behavior, health care availability and age.
Under the described conditions, the model dynamics can be written as: As natural births and deaths are not taken into account, the total number of individuals in the populations is considered to be constant, and the dynamics of the dead compartment can be written as a function of the other compartments: Besides, the recovered compartment is dependent on the susceptible and infected populations: 3) Therefore, the system dynamics can be rewritten without the removed and the dead variables, without any loss of generality: 3 System dynamics analysis

Equilibrium conditions
For the proposed SIRD model, equilibrium points can be obtained by applying the conditionsṠ = 0 andİ = 0. Then, it can be concluded that every point of the form (S * , 0) is an equilibrium point.
Equilibria at I = 0 are called disease-free, as opposed to equilibria at I = 0, which are called endemic [49]. The SIRD model presents the former but not the latter.

Equilibrium points classification
In order to analyze the local stability of these points, the Hartman-Grobman theorem is applied [50], and the general Jacobian (J) constructed as: (3.5) For (S, I) = (S * , 0) the Jacobian is reduced to: Hence, the eigenvalue set for the equilibrium points can be found by: Which gives: as the eigenvalue set. Considering that one of the eigenvalues is equal to zero, it is worth noticing that the eigenvalue set obtained from the equilibrium points allows a local stability analysis, with three possibilities: -If βS * N − (γ + Ω) > 0 ⇒ these fixed points are unstable [50], so there will be growth of cases for a small perturbation.
-If βS * N − (γ + Ω) < 0 ⇒ the Center-Manifold theorem must be applied in order to classify the stability condition.
Consequently, to proceed with the equilibrium stability analysis, all points satisfying βS * N − (γ + Ω) < 0 and S * 0 are studied applying the Center-Manifold theorem.
Considering S * > 0, a new parameter τ > 0 is introduced in order to translate the system so that the point (S* = τ ,0) is the origin (0,0) of the new system. This is done because the Center Manifold Theorem only refers to neighborhoods of the origin of the corresponding analyzed system. Consequently, the rewritten system, with Z = S − S * is given by: Consequently, studying the system (2.4) around (τ ,0) is the same as studying the new system (3.9) around (0,0). Nevertheless, the rewritten system doesn't possess the same physical interpretation as the SIRD model. In addition, substitute parameters are introduced in order to remove non-essential parameters regarding the interpretation of the stability condition. Then: Additionally, it is necessary to transform (3.11) into the form (A.18), as shown in Appendix A, by changing the basis into a basis of eigenvectors. Consequently: (3.12) Therefore, (3.13) To find h(u), as Appendix A shows, it is necessary to consider that the coordinates of any point in the center manifold must satisfy the conditions Now, after substituting (3.13) (with v = h(u)) into (3.15), it can be seen that h(u) = 0 is a solution of the partial differential equation (3.15). Then, the vector field restricted to the center manifold is given bẏ and, consequently, (S * , 0) is stable for the same set.
To complete the stability analysis, the case with S * = N (γ + Ω) β is considered, i.e., the Center-Manifold theorem [50] is applied for: In order to assess stability, a simulation is presented in As it can be seen, neighborhoods close to the origin with initial values I * < 0 are unstable, as the trajectories point towards negative infinite. However the system is only defined for I * > 0, considering that there no negative infected. Although Z can be negative as the original system was translated and even a small disturbance would never produces a state with I * < 0. Consequently, it is a degenerated equilibrium [51].

Basic reproduction rate
The basic reproduction rate R 0 is defined as the average number of secondary infections produced when one infected individual is introduced into a host population, where almost everyone is susceptible (S ≈ N ) [13]. For example, if the R 0 for COVID-19 in a population is 5, then each new case is expected to produce 5 secondary infections, assuming everyone around is susceptible. 8 Dieguez/Batistela/Piqueira Mathematically, R 0 is a threshold for stability of a disease-free equilibrium and is related to the peak and final size of an epidemic [49].

Phase portrait for different values of R 0
Figures (3) and (4) show the phase portrait of the system for different values of R 0 , specifically chosen to portrait the only two possible outcomes.
In (3), it is shown that, for R 0 > 1, the number of infected initially grows, but, asṠ < 0, at some point S * N eventually becomes < 1 R 0 , which implieṡ I < 0, and the trajectory always ends up in the horizontal axis (I = 0). This result is known as the herd immunity value. The interception point depends on the initial conditions, therefore, there can be more or less individuals that will never become infected in the course of the epidemic (final number of susceptible at disease-free equilibrium).
In Figure (4), with R 0 < 1, the system points toward the horizontal axis much faster, with the number of infected decreasing steeply all time. This was to be expected asİ < 0 for R 0 < 1.

Case studies and numerical simulations
In this section, two case studies are presented, with parameters derived from Brazil and Uruguay. These two countries were chosen because of their different basic reproduction number, as a result of particular population dynamics concerning the contact rate β. For Brazil, R 0 is estimated to be close to 2.82 [48] and for Uruguay, R 0 close to 1.13 [52]. These values were obtained using early infections data, as shown in figures (5) and (6), and are believed to present a scenario where there is no response to the pandemic, both individually or systemically.   Both graphs display a good fitting for initial data, with subsequent detachment between model and data, suggesting dynamical changes due to the implementation and effectiveness of isolation policies or individual confinement. In section 4.1 these set of parameters is used to simulate longer-term scenarios with no isolation that will serve as a reference for the evaluation of the control algorithm.

Simulation for Brazil
First, the model was simulated considering the parameters for Brazil, i.e., N = 209.5 million inhabitants, and I 0 = 10 individuals are considered, with the results shown in Figure 7. The chosen parameters conduct the simulation to a full blown-out pandemic spread. The last case ended about 120 days after the disease had begun being disseminated as we can see in figure 5.B. Figure 5.D shows a high mortality figure at the pandemic end, i.e., about 2.71% of the total population before the pandemic.   Figure 8 shows clearly the large proportion of the population that became infected at some point. At the end, less than 10% never became infected.

Simulation for Uruguay
The model was also simulated with parameters from Uruguay, i.e., N = 3.45 million inhabitants, considering I 0 = 10 individuals, and results shown in Figure 9.
As R 0 for Uruguay is lower than Brazil's, there is significant differences concerning the dissemination of the virus. Figure 9 shows that the pandemic lasted much longer in Uruguay than in Brazil, however, Uruguay had fewer dead individuals at the end, only about 0.64% of the total population before the pandemic. The low amount of individuals that ever got infected in Uruguay is highlighted in Figure 10, about 24% of the total population before the pandemic. Additionally, the lower number of infected suggests that having a lower R 0 due to lower β (contact rate) is a good indicator on what proportion the virus spread will achieve. This fact inspires the control algorithm used as reference of public policies employment, and is studied in the next section.

Spread control
Applications of Control Theory to Epidemiology range from Optimal Control, which can be used for social distancing [54], vaccine deployment policy [55], or even a mix of isolation and a vaccination program [56]; Model Predictive Control (MPC), which is also used to develop of social distancing [57,58] and vaccination policies [59]. However, the approach to be followed here is simpler, and aims to build an algorithm explicitly programmed through a feedback of state variables, similar to other heuristics [60,?] already studied. Practical usage of the control algorithm proposed should consider the availability and uncertainty of current infections data, which are usually embedded with underreporting and significant delays when compared to reality, specially in countries which do not have a thorough testing program organized. It should review the simplifications assumed in the model formulation in section 2 as well, addressing the peculiarities of the disease being analyzed.
The present algorithm should address, in the future, any disease that follows the same pattern of transmission of the coronavirus, that is, infections resulted from the transportation of contaminated droplets of secretions from infected individuals to susceptible ones and that suggests isolation is effective. The main difference will be the difficulty imposed by a new virus concerning the imposition of lower values of R t , parameter which will be described in next section. Viruses that are transmitted more easily impose a bigger challenge, as fewer contacts are needed to effectively pass the disease ahead.

Algorithm Implementation
Here, the implementation of a three-level stratified controller is studied, with the intention of mitigating the disease spread, aiming to maintain the demand of intensive care units (ICUs) lower than the available at all time, until a vaccine is available.
As it was shown, R 0 measures how the disease propagates when there are almost no infected individuals present (S ≈ N ). While R 0 provides valuable information on the viral dissemination dynamics when there is no developed immunity and information of the epidemic available, other factors begin to influence its dynamics during the spread course.
Therefore, the effective reproduction number R t , which describes the average value of secondary infections as a time-dependent function, is a more appropriate parameter to be analyzed amidst an epidemic outbreak [62]. In a homogeneous population, it is simple to obtain R t by multiplying R 0 and the susceptible proportion of the population at some instant t [63].
The idea to be developed is to regulate R t in order to maintain the health care resources at an acceptable level. The underlying mechanism is the manipulation of the contact rate of individuals β by cycling the imposition and relaxation of social distancing measures [46] as well as awareness campaigns to stimulate adoption of individual protection procedures.
Therefore β is no longer considered fixed, but time-dependent β(t). On the contrary, γ and Ω are still considered to be constant, as treatments such as the use of remdesivir or monoclonal antibodies, both approved by the American Federal Drug Administration (FDA), are prescribed only to severe cases whereupon the patient is already hospitalized and therefore, isolated [64,?] in the remdesivir case, while, for mono clonal antibodies, it is only approved for emergency use [66].
Levels of the controller are decided based on the current number of patients needing an ICU. This number is defined as 5% of the current number of infected individuals, assumed as the probability of someone who becomes infected to need intensive care [67].
Each level is composed of an interval from which R t can assume randomly any value. The interval is constructed with limits equal to ± 10% of the measured mean value. This is done to mimic the difficulty of imposing and tracking a precise value for R t for diseases with only a few epidemic studies. The stratified controller model was chosen because it is easier to implement than other forms of modeling social distancing, generally based on continuous values for R t , which are very hard to track and impose precisely. The algorithm is also adaptable structurally speaking as the R t should necessarily guarantee the desired behavior according to the model. The biggest difficulties rest in creating the appropriate public policies to reach the desired R t level.
Estimations of the economic impact of COVID-19 pandemic situate close to trillion dollars figures globally [73,74,75], with devastating shocks on various industries, and the imposition of social distancing policies should have been, theoretically, a big part of this loss. However, there is significant evidence that quarantine is effective against not implementing containment mechanisms [76], more effective when applied early [77], and that deterioration of economic conditions preceded the introduction of isolation policies [78] or that the culprit of COVID-19 recession is COVID-19 itself [79]. Therefore, practical usage of the three-level controller should lead to economic benefits, not to mention saving lives should be worth the cost nevertheless.
The three-level controller is described in Table II, with the calibration level (mitigation or suppression) related to the value of reproduction number associated. Mitigation and suppression measures differ in whether they aim to reduce the reproduction number, R t , to less than 1 (suppression) or to merely slow spread by reducing R t , but not to less than 1 (mitigation) [80]. Therefore, Level A can be described as mitigation strategy and Level C as a suppression one, generally called a lock down, while Level B situates between both.

E[Rt]
Measures taken Value source When?
Self-imposed measures are stimulated in order to accelerate awareness spread. Prevention measures such as mask-wearing, handwashing and self-imposed social distancing can be described as reductions in infectious output, susceptibility and contact rate, respectively [81].
As such measures stack up additively [81], a 10% efficacy on each one reduces the effective contact rate by 30%.
Intensive care unit demand is lower than 10% of the total amount available. All the actions taken at a certain point will last for a regular chosen duration, during all the course of the pandemic. This is done in order to reduce social and economic uncertainty, as the population have time to comply.

Control strategy results
In this section, the three-level control strategy, described in the former section is simulated assuming an hypothetical scenario for an isolated city in Brazil with 10 infected initial cases, a population of 100,000 inhabitants and R 0 the same as the whole country (2.82). This number was maintained for the 15 initial days to replicate initial unawareness. Brazil is chosen due to its high baseline contact rate, represented by the high value of R 0 , as opposed to Uruguay. Therefore, as the number of cases and individuals in need for ICUs grow faster than the system can adapt, it represents a need for the establishment of isolation policies.
The number of ICUs available is considered to be 50, very close to the proportion of beds/individuals obtained for the state of Sao Paulo before the pandemic, in 2018 [68].
The total period analyzed is one year, because this is the period assumed for the development of a vaccine, with reference to the Pfizer vaccine, which started being developed January 2020 with the release of the SARS-CoV-2 genome [69], and initially applied in the UK at the beginning of December 2020 [70].
Three different strategies were considered: no attempt to control the disease spread; updating R(t) every 30 days; and updating R(t) every 21 days.

No attempt to control the spread of the disease
As expected, figure 11 shows the potential of COVID-19 to quickly overwhelm the health care system, with ICU demand surpassing the availability of beds less than two months after the beginning of the pandemic. The randomness added did not result in significant differences with respect to the baseline established in section 4. number R t for each cycle. It can be seen that, after periods of Level C (lock down), ICU demand went below 10% of the total number of intensive care beds available, which made the algorithm choose Level A for the following cycle.
Every time that this happened, the control failed to prevent health care collapse, as 30 days proved to be enough time for Level A to allow a large increase of demand for ICU, due to the exponential growth nature of the spread.  The improvement over the baseline was very small, as can be seen in figure  12.D, with then number of deceased at the end of the period being very close to the number of deceased at the end in the baseline.
Therefore, in order to prevent the ICU demand from surpassing the supply of intensive care beds, the period between updates of the algorithm is reduced, and results shown in the next sub-section.

Updating R(t) every 21 days
For a 21 days interval between updates, the main goal of guaranteeing appropriate care for COVID-19 patients is attained. During the 360 days of simulation, there was no instant at which the demand for intensive care unit was bigger than the supply as shown by figure 15, and only three periods of Level C (lock down), allowing the individuals to have greater freedom without sacrificing effectiveness. Figure 16 supports the effectiveness of the algorithm with 21 days cycle, showing much fewer deaths at the end of the simulation than in simulations I (baseline) and II (30 days cycle).

Conclusions
Inspired by the different policies of social distances adopted by Brazil and Uruguay, showing remarkable differences concerning the results on controlling COVID-19 epidemic, a framework for modeling and applying public policies in the form of social distancing was proposed aiming the reduction of the disease spread.
A three-level controller was proposed and simulations indicate that a 21day update strategy has shown good results, avoiding the health care system to collapse and presenting much fewer deaths at the end of the process.
To improve the model, it could be considered a varying R t for each level and different number of levels. Uncertainties and delays could improve the observability of the number of infected cases.
Despite being an intuitively obvious conclusion, it was shown that not trying to mitigate the spread of the disease is not a wise option. The algorithm with 21 days between updates presents, for the example, almost 650 deaths at the end of the year, whereas doing nothing resulted in approximately 2700.
Therefore, it is possible to say that the framework outlined is a good and simple reference model to be followed when designing techniques to address the COVID-19 disease spread, and possibly other diseases that follow the same pattern of transmission as the corona virus.