A quantitative framework for exploring exit strategies from the COVID-19 lockdown
AA quantitative framework for exploring exit strategies from the COVID-19 lockdown
A.S. Fokas a,b , J. Cuevas-Maraver c,d and P.G. Kevrekidis e,f a Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, U.K. b Department of Civil and Environment Engineering, University of Southern California, 90089, Los Angeles, Ca, USA c Grupo de Física No Lineal, Departamento de Física Aplicada I, Universidad de Sevilla. Escuela Politécnica Superior, C/ Virgen de África, 7,41011-Sevilla, Spain d Instituto de Matemáticas de la Universidad de Sevilla (IMUS). Edificio Celestino Mutis. Avda. Reina Mercedes s/n, 41012-Sevilla, Spain e Department of Mathematics and Statistics, University of Massachusetts, Amherst, MA 01003-4515, USA f Mathematical Institute, University of Oxford, Oxford, UK
A R T I C L E I N F O
Keywords :COVID-19 ModelingOrdinary Differential EquationsCumulative Death ModelingParameter Estimation and Optimiza-tionLockdown Exit Strategies
A B S T R A C T
Following the highly restrictive measures adopted by many countries for combating the currentpandemic, the number of individuals infected by SARS-CoV-2 and the associated number ofdeaths steadily decreased. This fact, together with the impossibility of maintaining the lockdownindefinitely, raises the crucial question of whether it is possible to design an exit strategy basedon quantitative analysis. Guided by rigorous mathematical results, we show that this is indeedpossible: we present a robust numerical algorithm which can compute the cumulative number ofdeaths that will occur as a result of increasing the number of contacts by a given multiple, usingas input only the most reliable of all data available during the lockdown, namely the cumulativenumber of deaths.
1. Introduction
On December 31, 2019, the Chinese government reported a cluster of pneumonia cases of unknown cause that waslater identified as a result of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). This is one of themost serious manifestations of an infection and associated disease (termed COVID-19) caused by a coronavirus. Likeearlier outbreaks caused by two other pathogenic human respiratory coronaviruses, namely, the severe respiratorysyndrome coronavirus (SARS-CoV) and the Middle East respiratory syndrome coronavirus (MERS-CoV), SARS-CoV-2 causes a respiratory disease that it is often severe. In addition, SARS-CoV-2 can attack many vital organs ofthe body, and can also lead to severe neurological disorders, including the Guillain-Barré syndrome [1, 2]. Fortunately,SARS-CoV-2 is associated with lower mortality than its above predecessors; however, it is more contagious [3]. As aresult of this fact and the lack of early measures for curtailing its spread, it has caused a pandemic, which is consideredas the most serious threat to public health since the pandemic caused by the 1918 influenza (the 1918 pandemic, whichis the deadliest event in human history, caused more than 50 million deaths which corresponds to 200 million deathsin todayâĂŹs population).The scientific community is playing a crucial role in combating the above threat: from elucidating fundamentalfeatures of SARS-CoV-2 and mechanisms of its transmissibility [4], to addressing the vital question of a pharmaco-logical treatment and the development of an effective vaccine [5]. For example, the viral genome of SARS-CoV-2has been sequenced [6]. Also, mechanisms underlying the increased transmissivity of SARS-CoV-2 have been tracedto its dual receptor attachment in the host cells. In particular, it has been shown that the attachment of the virus tothe surface of respiratory cells is mediated by certain viral proteins which bind not only to the angiotensin convert-ing enzyme-2 (ACE-2) receptor [7], but also to sialic acid containing glycoproteins and gangliosides that reside oncell surfaces [8] . Regarding pharmacological interventions, a randomized, controlled study involving hospitalizedadult patients with confirmed SARS-CoV-2 infection, showed that there was no benefit from the antiviral regime oflopinavir-ritovir (which is an effective treatment in patients infected with human immunodeficient virus) [9]. Sim-ilarly, the combination of the anti-malarial medication chloroquine, and its derivative hydroxychloroquine, with orwithout the antibiotic azithromycin, may not only be ineffective, but also potentially harmful [10, 11, 12] (see also themeta-analysis [13]). ORCID (s): This is to be contrasted with the case of SARS-CoV, which binds only to ACE-2 receptors [7, 8].
A.S. Fokas, J. Cuevas-Maraver and P.G. Kevrekidis:
Preprint submitted to Elsevier
Page 1 of 14 a r X i v : . [ q - b i o . P E ] A ug quantitative framework for exploring exit strategies from the COVID-19 lockdown In addition, the scientific community has had a vital impact in decisions made by policy makers of possible non-pharmacological approaches to limit the catastrophic impact of the pandemic. For example, two possible strategies,called mitigation and suppression, are thoroughly discussed in [14, 15]; in the early stages of the pandemic, UKwas following mitigation, but after the publication of this report, it adopted suppression. This policy, which wasalready implemented in several other countries (of course there were also notable exceptions such as Sweden and Brazilwith serious associated consequences), has led to the curtailing of the pandemic. Indeed, the number of individualsreported to be infected with SARS-CoV-2, as well as the number of related deaths, steadily decreased in the countriesthat adopted severely restrictive measures, known as âĂŸlockdownâĂŹ; see e.g. the data in [16]. This welcomedevelopment, together with the impossibility of maintaining indefinitely the lockdown conditions (both for obviouseconomical reasons and for the induced serious psychological effects) raises the following question: is it possible todesign an exit strategy based on the quantitative analysis of the effect of easing the lockdown measures? The purposeof this work is to provide an affirmative answer to this question, concentrating on the number of associated deaths, andthe unique identification for a canonical epidemiological model of parametric combinations that determine it.We note in passing that the use of mathematical methods in a plethora of biological problems at the level of (bothdiscrete and continuum) model analysis and at that of finding a wide variety of special solutions is a theme of intenseinterest in recent years. As some relevant examples, we mention the works of [17, 18] (while similar methods have beenused in nonlinear engineering and mathematical physics problems [19, 20, 21, 22]). In the recently very highly activefront of COVID-19 modeling, some of the efforts have been directed at modeling the early stages of the pandemic [23];others have focused on designing a pandemic response index (to quantify/rank the response of different countries) [24]or towards quantifying the response of different regions within a country, e.g., the states within the USA [25]. Similarly,models have focused on cruise ships [26], on cities [27], as well as states/provinces [27, 28, 29, 30], but also variouscountries [15, 31, 32, 33, 34, 35], aside from the prototypical examples of Wuhan, China [36], and some among thehard-hit Italian provinces [37].Before providing details of our proposed methodology for the computation of the number of deaths following aspecific increase of contacts between asymptomatic individuals infected with COVID-19, it should be emphasizedthat the answer to the above question is literally a matter of life and death: (i) No therapeutic intervention has beenproven so far effective for the treatment of the severe illnesses and side-effects caused by SARS-CoV-2. (ii) Reportedmortality rates differ drastically between different countries and are crucially affected by age. For example, in the largeststudy from China involving 1099 hospitalized patients with laboratory confirmed SARS-CoV-2 infection of medianage 47 years, only 5% needed admission to the intensive care unit, 2.3% underwent invasive mechanical ventilation,and 1.4% died [38]. On the other hand, following the identification on February 28, 2020 of a confirmed case ofCOVID-19 in a nursing facility in Washington, USA, as of March 18, 101 residents of this facility and 50 healthcare personnel were confirmed with COVID-19, and in addition 16 infected visitors were epidemiologically linkedwith this facility. Hospitalization rates for residents, staff, and visitors, were 54.5%, 6%, and 50% respectively; thecorresponding mortality rates were 33.7%, 0%, and 6.2% [39]. (iii) Despite the fact that the unprecedented efforts forthe development of a vaccine take place within a framework of explosive progress in basic scientific understanding thathas occurred in the areas of genomics and structural biology, it is not expected that a vaccine will be available for atleast several months [40]. From the above remarks it becomes clear that the lockdown measures must be eased withoutthe benefit of any substantial pharmacological cover, which is desperately needed especially for older persons (andindividuals with a variety of diseases such as Diabetes type 2, hypertension, and respiratory disorders), and withoutthe anticipation of the imminent availability of a vaccine.There is no doubt that the solution to the vitally important problem of âĂŸhow to defeat the current pandemicâĂŹ,will finally be provided by medicine and biology via the development of appropriate pharmacological treatments andan effective vaccine. However, it appears that at the moment, there exists a unique opportunity for mathematicalmodeling and associated analysis to contribute definitively to the partial addressing of this problem. In particular,SIR (susceptible, infected, recovered) type models are widely accepted in mathematical epidemiology [41]; and thistype of models can be modified to capture some key features of the present pandemic (such as the crucial role ofasymptomatically infected individuals [26, 42, 43]). Thus, it follows that the rigorous analysis of such a class ofmodels provides a possible approach towards studying quantitatively the effect of easing the lockdown measures. Ifsuch analysis can be leveraged to give rise to accurate and reproducible computational results based on reliable data,then it would be possible to explore systematically the design of a safe lockdown exit.The presentation of our efforts in the above direction is structured as follows. In section 2, we provide an overviewof the computational algorithm that we propose, as well as its epidemiological implications towards identifying the
A.S. Fokas, J. Cuevas-Maraver and P.G. Kevrekidis:
Preprint submitted to Elsevier
Page 2 of 14 quantitative framework for exploring exit strategies from the COVID-19 lockdown model parameters and also evolving the model in the future towards release of lockdown conditions. Then, in section3, we delve into the details of our mathematical methods and corresponding computational results. Finally, in section4, we summarize our findings and discuss a number of future recommendations and interesting research directions.
2. A computational algorithm based on a rigorous mathematical result
The model adopted in this work, which is discussed in detail in the Methods and Results section, involves 6 (firstorder) ordinary differential equations (ODEs) uniquely specified by 9 constant parameters (and, of course, well-posedunder the introduction of associated initial conditions). One of the relevant parameters, namely the constant 𝑐 ,denotes the effect of the interaction of the asymptomatic, infected individuals with those susceptible to be infected. Thisparameter is particularly important for our analysis: the easing of the lockdown conditions would result in increasingthe number of contacts (to which 𝑐 is proportional) among asymptomatic and susceptible individuals, and this effectcan be straightforwardly incorporated in the model by replacing 𝑐 with 𝜁𝑐 , where 𝜁 is a fixed number, such as or ; this will be referred to, respectively, as the doubling or tripling of the number of contacts.How can our model yield a computational approach to designing a safe lockdown exit? Suppose that, somehow,given appropriate data obtained during the lockdown period, the values of the parameters specifying the model couldbe determined (and the associated initial conditions prescribed). Then, using these values, replacing 𝑐 by 𝜁𝑐 , andsolving the resulting 6 ODEs, the number of infected, hospitalized, recovered, and deceased individuals in the post-lockdown period could be computed. It turns out that the above scenario, which would yield information about allthe basic features of the post-lockdown state, is impossible. An implicit assumption in this class of epidemiologicalmodels is that the model parameters can be uniquely determined from an appropriate set of data. However, this isfar from a trivial assumption and the whole branch of identifiability of the models is associated with this issue [44].In this connection it should be noted that the problem of determining the model-parameters from the given data canbe formulated as an inverse problem, and such problems are notoriously difficult, especially regarding the importantquestion of uniqueness; namely, proving that the given data give rise to a unique set of parameters . Actually, we haveshown that it is impossible to determine uniquely all 9 model parameters from a reliable set of given data. This is indeedan example of dimension reduction: this represents the crucial aspect of whether a given model outcome depends noton individual parameters alone (except for one of the model parameters) but rather on combinations formed by suitable(irreducible) combinations of parameters; for a recent, data-driven example, see [46]. In addition to the existence ofthe above prohibitive result, many of the needed data are unavailable. For example, although the model involves thetotal number of infected individuals, the available data are for the âĂŸreported infectedâĂŹ individuals.By concentrating on the number of deaths, we have been able to overcome both of the above difficulties. Indeed,we have shown that it is possible to determine the single model parameter and the 6 specific combinations of the 9model parameters that characterize a certain 4th order ODE specifying the time-evolution of the number of deaths.Furthermore, we have established that this can be achieved uniquely by using the most reliable of all available data,namely, the data for the cumulative number of deaths.The above mathematical results give rise to the following algorithm (which we have showcased in a number of selectexamples): starting with the death data during the lockdown period, we compute uniquely the single model parameteras well as the 6 combinations of the original 9 model parameters that specify the 4th order ODE determining theevolution of the number of deaths. The parameter 𝑐 enters one of these combinations, and fortunately, it enters in ahomogeneous manner; thus, replacing 𝑐 by 𝜁𝑐 results in replacing the relevant combination by 𝜁 times its originalvalue. Then, using this combination, the single model parameter, and the remaining 5 combinations determined fromthe death data, we can proceed to forward time-step the resulting 4th order ODE. This uniquely yields the number ofdeaths in the post-lockdown period (the concrete procedure used to test the robustness of our algorithm, is describedin the Methods and Results section).We applied the above approach to the epidemics of the countries of Portugal and Greece, as well as the autonomouscommunity of Andalusia in Spain. These have been selected due to their notable similarities (geographic location, aswell as similar populations in the vicinity of 8-10 million residents), but also due to their significant differences inconnection to the response to the pandemic. Greece has had an extremely low tally of deceased individuals as a resultof the pandemic, being one of the notable success examples of early application of lockdown measures [47]. Portugal For example, it is rigorously established in [45] that different neuronal electric currents give rise to the same data obtained via electroen-cephalographic (EEG) recordings. Thus, it is impossible to obtain uniquely the current via the solution of the inverse problem associated withEEG.
A.S. Fokas, J. Cuevas-Maraver and P.G. Kevrekidis:
Preprint submitted to Elsevier
Page 3 of 14 quantitative framework for exploring exit strategies from the COVID-19 lockdown is also a case with relatively low cumulative numbers, although, at around the same population as Greece, it currentlyhas more than 7 times the number of deaths. Andalusia, on the other hand, has a smaller population (by about 2 million)than the other 2, but has already reported deaths (i.e., more than Portugal and more than 8 times more thanGreece). Hence, these are interestingly distinct examples. In the case of Greece, and through the findings reflecting thelow transmissivity of the virus, we find that a doubling of the contacts in the entire population, i.e., 𝜁 = 2 , would onlyslightly increase the number of deaths from 158 to 167. On the other hand, the effect of such a widespread easing of thelockdown measures would be quite different in the case of Portugal (from 1433 to 40727) and of Andalusia (from 1451to 26846). Furthermore, far more dramatic is the impact in the case of changing 𝜁 to . In that case, the number ofdeaths is drastically increased even in Greece (18474); the situations in Portugal (84014) and Andalusia (106855) maybecome catastrophic. The derivation of the above results is the consequence of a stable numerical algorithm capturingthe key finding presented here, namely, that it is possible to: (i) identify the irreducible sets of parameters associatedwith the system; and (ii) to subsequently utilize (only) those combinations towards the forward prediction of the deathtally upon different scenarios of easing the lockdown measures.
3. Methods and Results
Let 𝐼 ( 𝑡 ) denote the infected (but not infectious) population. An individual in this population, after a median 5-dayperiod (required for incubation [49]) will either become sick or will be asymptomatic, as shown in the flowchart ofFig. 1 , which represents all the relevant populations and the transitions between them. The sick and asymptomaticpopulations will be denoted, respectively, by 𝑆 ( 𝑡 ) and 𝐴 ( 𝑡 ) . The rate at which an infected person becomes asymptomaticis denoted by 𝑎 ; this means that each day 𝑎𝐼 ( 𝑡 ) persons leave the infected population and enter the asymptomaticpopulation. Similarly, each day 𝑠𝐼 ( 𝑡 ) leave the infected population and enter the sick population. The asymptomaticindividuals recover with a rate 𝑟 , (i.e., similarly each day 𝑟 𝐴 ( 𝑡 ) leave the asymptomatic population and enter therecovered population), which is denoted by 𝑅 ( 𝑡 ) . The sick individuals either recover with a rate 𝑟 or they becomehospitalized 𝐻 ( 𝑡 ) with a rate ℎ . In turn, the hospitalized patients also have two possible outcomes of the medicalintervention efforts; either they recover with a rate 𝑟 , or they become deceased, 𝐷 ( 𝑡 ) , with a rate 𝑑 . I ! " a AS H RD s ! h d ! $ % $ A% S Figure 1:
Flowchart of the populations considered in the model and the rates of transformation between them. Thecorresponding dynamical equations are Eqs. (1)–(6).
It is straightforward to write the above statements in a mathematical form; this gives rise to the equations below: 𝑑𝐴𝑑𝑡 = 𝑎𝐼 − 𝑟 𝐴, (1) 𝑑𝑆𝑑𝑡 = 𝑠𝐼 − ( ℎ + 𝑟 ) 𝑆, (2) 𝑑𝐻𝑑𝑡 = ℎ𝑆 − ( 𝑟 + 𝑑 ) 𝐻, (3) 𝑑𝑅𝑑𝑡 = 𝑟 𝐴 + 𝑟 𝑆 + 𝑟 𝐻, (4) 𝑑𝐷𝑑𝑡 = 𝑑𝐻. (5)In order to complete these equations, it is necessary to describe the mechanism via which a person can becomeinfected. For this purpose, we follow the standard assumptions made in the SIR-type [41] epidemiological models: let An interval of 3-10 days captures 98% of the cases.
A.S. Fokas, J. Cuevas-Maraver and P.G. Kevrekidis:
Preprint submitted to Elsevier
Page 4 of 14 quantitative framework for exploring exit strategies from the COVID-19 lockdown 𝑇 denote the total population and let 𝑐 be the transmission rate proportional to the number of contacts per day made byan individual with the capacity to infect. Such a person belongs to 𝑆 , or 𝐴 , or 𝐻 . However, for simplicity we assumethat the hospitalized population cannot infect; this assumption is based on two considerations: first, the strict protectivemeasures taken in the hospital, and second, the fact that hospitalized patients are infectious only for part of their stayin the hospital. The asymptomatic individuals are (more) free to interact with others, whereas the (self-isolating) sickpersons are not; yet, the viral loads of the two have been argued to be similar [42]. Thus, we use 𝑐 and 𝑐 to denote thecorresponding transmission rates (per person in the respective pools) of the asymptomatic and sick respectively. Thenumber of persons available to be infected is 𝑇 − ( 𝐼 + 𝑆 + 𝐴 + 𝐻 + 𝑅 + 𝐷 ) . Indeed, the susceptible individuals consistof the total population minus all the individuals that after going through the course of some phase of infection, eitherbear the infection at present ( 𝐼 + 𝐴 + 𝑆 + 𝐻 ) or have died from COVID-19 ( 𝐷 ) or are assumed to have developed(even if temporary, but sufficient for our time scales of interest) immunity to COVID-19 due to recovery ( 𝑅 ). The rateby which each day individuals enter 𝐼 is given by the product of the above expression with 𝑐 𝐴 + 𝑐 𝑆 . At the sametime, as discussed earlier, every day ( 𝑎 + 𝑠 ) 𝐼 persons leave the infected population. Thus, the rate of change of 𝐼 ( 𝑡 ) reads: 𝑑𝐼𝑑𝑡 = [ 𝑇 − ( 𝐼 + 𝑆 + 𝐴 + 𝐻 + 𝑅 + 𝐷 )] ( 𝑐 𝐴 + 𝑐 𝑆 ) − ( 𝑎 + 𝑠 ) 𝐼. (6)The above model depends on the given (total population) constant 𝑇 and on the 9 parameters 𝑑 , ℎ , 𝑠 , 𝑎 , 𝑐 , 𝑐 , 𝑟 , 𝑟 , 𝑟 . The fundamental inverse problem that we consider herein is the following: which specific parameters and combi-nations of the model parameters can be uniquely determined from the knowledge of the death data? We utilize herethe data on the deceased individuals because it is the most reliable time series among the ones available. Indeed,there exists a significant uncertainty regarding the number of true infections (vs. the ones officially reported). Thisuncertainty is even more significant with respect to the current situation regarding the asymptomatics which, whileconcretely studied in some cases [26], presents tremendous variability [43] between studies. As indicated above, theissue at hand is crucial from the point of view of modeling, both regarding issues of dimension reduction [46], as wellas those of identifiability of the models [44].It turns out that the cumulative number of deceased, 𝐷 , satisfies a 4th order nonlinear ODE, which is uniquelydetermined by the constants 𝛼 = ℎ𝑠𝑑 ( 𝑇 − 𝜇 ) and 𝛽 , where 𝜇 and 𝛽 are integration constants, the model parameter 𝑟 , as well as the following 5 combinations of the model-parameters: 𝐶 = 𝑎𝑐 ∕( ℎ𝑠𝑑 ) , 𝐶 = 𝑐 ∕( ℎ𝑑 ) , 𝐹 = 𝑎 + 𝑠 , 𝑅 = 𝑟 + ℎ , 𝑅 = 𝑟 + 𝑑 . It is relevant to note here, in passing, that the initial condition dependent 𝜇 ≪ 𝑇 practically,and hence it is possible through 𝛼 , in principle, to provide a good approximation of the product ℎ𝑠𝑑 . This basic ODEassumes the form 𝑞 + 𝑞 ( 𝑞 + 𝑟 ln( | 𝑞 | ) ) = 0 , (7)where 𝑞 = ( 𝐶 + 𝐶 ) 𝐷 (2) + [ 𝐶 ( 𝑅 + 𝑅 ) + 𝐶 ( 𝑅 + 𝑟 ) ] 𝐷 (1) + 𝑅 ( 𝐶 𝑅 + 𝑟 𝐶 ) 𝐷 + 𝛽,𝑞 = 𝐷 (3) + 𝑘 𝐷 (2) + 𝑘 𝐷 (1) + 𝑘 𝐷 − 𝛼,𝑞 = 𝐷 (4) + 𝑘 𝐷 (3) + 𝑘 𝐷 (2) + 𝑘 𝐷 (1) , (8)with superscripts denoting the number of derivatives with respect to 𝑡 , and the constants appearing in Eq. (8) definedas follows: 𝐶 = 𝑎𝑐 ℎ𝑠𝑑 , 𝐶 = 𝑐 ℎ𝑑 , 𝑅 = 𝑟 + ℎ, 𝑅 = 𝑟 + 𝑑, (9) 𝑘 = 𝐹 𝑅 𝑅 , 𝑘 = 𝐹 ( 𝑅 + 𝑅 ) + 𝑅 𝑅 , 𝑘 = 𝐹 + 𝑅 + 𝑅 , 𝐹 = 𝑎 + 𝑠. We next explain how the central equation (7) has been obtained and what the practical implications of the aboveresults are before turning to the examination of our numerical findings. The main idea is to work backwards andutilize Eq. (5) as a way to express the 𝐻 ( 𝑡 ) in terms of 𝐷 (1) , then Eq. (3) to express 𝑆 ( 𝑡 ) in terms of 𝐷 (1) and 𝐷 (2) ,then Eq. (2) to express 𝐼 ( 𝑡 ) in terms of 𝐷 (1) , 𝐷 (2) and 𝐷 (3) ; do the same for 𝐴 ( 𝑡 ) + 𝑅 ( 𝑡 ) etc. Reaching back at thelevel of Eq. (6), one possible direction is to solve this equation as an algebraic one for 𝐴 and then substitute back toEq. (1). This leads to an interesting in its own right nonlinear 5th order ODE. However, one can obtain an even moresubstantial and quite remarkable simplification, as a result of the structure of the system. It turns out that the quantity A.S. Fokas, J. Cuevas-Maraver and P.G. Kevrekidis:
Preprint submitted to Elsevier
Page 5 of 14 quantitative framework for exploring exit strategies from the COVID-19 lockdown [ 𝑇 − ( 𝐼 + 𝑆 + 𝐴 + 𝐻 + 𝑅 + 𝐷 )] can be expressed in terms of the 𝐷 ( 𝑡 ) time series according to the form provided by 𝑞 above; it is a direct calculation to see that 𝑑𝐼 ∕ 𝑑𝑡 + ( 𝑎 + 𝑠 ) 𝐼 is essentially the derivative of this expression that entersinto 𝑞 . As a result, the ratio of these expressions 𝑞 ∕ 𝑞 becomes a logarithmic derivative, enabling us, when utilizedtogether with Eq. (1), to identify a âĂŸfirst integralâĂŹ of the relevant equation that ultimately results into Eq. (7) at themere expense of introducing an additional integration constant ( 𝛽 ). In fact, in addition to the remarkable simplificationof obtaining this as a 4th order ODE system, the presence of 𝛽 offers a useful benchmark for the numerical method aswe will see below.There is an additional, important mathematical observation that concerns the stability calculation of the healthystate of the system which in our model is represented by ( 𝐼, 𝑆, 𝐴, 𝐻, 𝑅, 𝐷 ) = (0 , , , , , i.e., the null state. Thespectral stability analysis of this state determines whether an epidemic will grow or decay on the basis of the dominantstability eigenvalue 𝜆 of the relevant state. The destabilization threshold 𝜆 = 0 corresponds to a basic reproductionnumber 𝑅 = 1 [41]. Positive 𝜆 leads to 𝑅 > and spreading of the epidemic, whereas 𝜆 < leads to the epidemicsubsiding over time ( 𝑅 < ). Starting from a healthy state in which case 𝜇 = 0 and 𝛼 = ℎ𝑠𝑑𝑇 in the expressionsabove, we find that the relevant eigenvalue can be obtained by the largest eigenvalue of the polynomial: 𝜆 + 𝜆 ( 𝐹 + 𝑟 + 𝑅 ) + 𝜆 [ 𝐹 ( 𝑟 + 𝑅 ) + 𝑟 𝑅 − 𝛼 ( 𝐶 + 𝐶 ) ] + 𝐹 𝑟 𝑅 − 𝛼 ( 𝐶 𝑅 + 𝐶 𝑟 ) = 0 . (10)Remarkably, the relevant dominant eigenvalue depends solely on the reduced quantities (including, of course, the roleof the total population introduced through 𝛼 ) of the original system. In fact, the direct application of the Routh-Hurwitzconditions [51] leads to the criterion 𝛼 ( 𝐶 𝑅 + 𝐶 𝑟 ) < 𝐹 𝑟 𝑅 , (11)for the subsiding of the epidemic (while the opposite sign of the inequality will lead to its spreading, and equalitywill lead to a vanishing eigenvalue or 𝑅 = 1 ). The above result shows that 𝑟 and the reduced set of parametriccombinations identified earlier, not only determine the most important tally of the impact of the epidemic, namely thenumber of deaths, but also the potential capacity of the epidemic to spread.It is important to mention that a similar analysis regarding the evolution of 𝐷 ( 𝑡 ) can be performed in the case oftwo subpopulations, young and older [48]. Then, Eq. (7) is replaced by the two equations 𝑞 𝑙 + 𝑃 𝑞 𝑙 = 0 , (12)where 𝑞 𝑙 and 𝑞 𝑙 , are defined via the equations obtained from (8) by simply inserting superscripts y or o, whereas P isdefined by: 𝑃 = ∑ 𝑙 = 𝑦,𝑜 𝐶 𝑙 𝐴 𝑙 + 𝐶 𝑙 𝑆 𝑙 . (13)Here, 𝐴 𝑙 and 𝑆 𝑙 correspond to the asymptomatic and sick populations of each pool (young or older) that can bealgebraically expressed in terms of 𝐷 ( 𝑡 ) of the two sub-populations and their derivatives.Our theoretical analysis has provided the lumped (reduced) sets of coefficients in the spirit of dimension reductionthat can be uniquely identified on the basis of the time series of 𝐷 ( 𝑡 ) , namely 𝑟 , 𝐶 = 𝑎𝑐 ∕( ℎ𝑠𝑑 ) , 𝐶 = 𝑐 ∕( ℎ𝑑 ) , 𝐹 = 𝑎 + 𝑠 , 𝑅 = 𝑟 + ℎ , 𝑅 = 𝑟 + 𝑑 . Additionally, one can identify the two integration constants 𝛼 and 𝛽 . In thepractical implementation of our algorithm we will utilize the knowledge of the existence of such an “optimal” set of(lumped) parameters to leverage the ODE of Eq. (7) to identify these parameters. The idea is that given a prescribedtime-series for 𝐷 ( 𝑡 ) and its derivatives, the left- hand side of Eq. (7) can be evaluated. For a perfect time-series andan exact integration procedure, the unique global optimum parameter set would render the time series of this left-handside (that we will refer to as 𝑓 hereafter) vanishing. As a measure of its vanishing, since practically the time seriesamounts to a vector at given times, we will use the || 𝑓 || ≡ || 𝑓 || 𝑙 ; other norms can be used as well.Despite the existence of the above well-defined mathematical procedure, there arise multiple practical complica-tions. In particular, even if one has a âĂŸperfectâĂŹ time series there is still the approximate nature of numericalintegration. In what follows, we will generate such a perfect series, that we will refer to as “synthetic” hereafter, viathe integration of Eqs. (1)-(6) for given (prescribed) parameters. In practice, computations have a finite numericalaccuracy that will be reflected in the approximate as opposed to the exact vanishing of || 𝑓 || . Naturally, the situationwith practical data-based time series, is far more difficult. In that case, we will only know 𝐷 and we need to infer its A.S. Fokas, J. Cuevas-Maraver and P.G. Kevrekidis:
Preprint submitted to Elsevier
Page 6 of 14 quantitative framework for exploring exit strategies from the COVID-19 lockdownParameter Exact Finite differences Sigmoid 𝐶 . −5 . −5 . −5 𝐶 . −6 . −6 . −6 𝑅 . . . 𝑅 . . . 𝐹 . . . 𝛼 . . . 𝑟 . . . 𝛽 −1 . . . Table 1
Lumped parameter combinations (first column) of the local, but also global within our considerations, minima for thesynthetic time series example and their numerical values for the exact up to numerical accuracy case (second column), thefinite difference derivative approximation (third column) and the sigmoid approximation (fourth column) approaches. 𝐶 𝐶 || 𝑓 || .
47 × 10 −9 . . .
13 × 10 −5 . . .
04 × 10 −5 . . .
16 × 10 −4 . . .
79 × 10 −4 Table 2
For the exact (up to numerical accuracy) case, the value of the global optimum and of the “next best things” (first twocolumns) and its corresponding fitness (third column) are given. The first two columns are given in terms of multiples ofthe ( 𝐶 , 𝐶 ) pair provided in Table 1. (noisy) first four time-derivatives, before substituting all the relevant time series ( 𝐷, 𝐷 (1) , 𝐷 (2) , 𝐷 (3) and 𝐷 (4) ) into 𝑓 ,in order to perform the minimization.We explored two possibilities in order to examine whether our approach could be brought to bear in practical situa-tions. In the first one, we considered finite difference-based approximations of the prescribed time series for 𝐷 ( 𝑡 ) ; 𝐷 ( 𝑡 ) constitutes the only provided reliable data to be used in the inverse-problem formulation for the identification of therelevant parameters. In the second approach we focused on a logistic formula approximation earlier explored system-atically, e.g., in [52, 53]. In the latter case, we can consider the variable 𝑥 = 𝐷 ∕ 𝐷 𝐹 (where 𝐷 𝐹 is the (asymptotically)final number of deceased individuals); then , the differential equation (7) becomes 𝑥 (1) = 𝑘𝑥 (1 − 𝑥 ) ⇒ 𝑥 (2) = 𝑘 ( 𝑥 − 3 𝑥 + 2 𝑥 ) ⇒ 𝑥 (3) = 𝑘 ( 𝑥 − 7 𝑥 + 12 𝑥 − 6 𝑥 ) ⇒ 𝑥 (4) = … (14)The advantage of the second methodology is its simplicity when employed in practical data-driven settings, as wellas the smoothness of the quantities involved. Namely, the given time series for 𝐷 ( 𝑡 ) is fitted into a sigmoid form: 𝐷 ( 𝑡 ) = 𝐷 𝐹 ∕(1 + 𝑏𝑒 − 𝑘𝑡 ) , where 𝑏 (for given 𝐷 𝐹 ) is essentially related to the initial condition 𝐷 (0) = 𝐷 𝐹 ∕(1 + 𝑏 ) and 𝑘 is the truly essential (and most robust [52, 53]) piece of information of the sigmoid fit that enters the differentialequations of (14). Substitution of these ODEs into Eq. (7) ultimately leads to a time series 𝑓 , or upon expansion (usingthe dominance of the term involving 𝑇 within the logarithmic term) to a polynomial. The minimization of this timeseries (in what is shown below) or of the associated polynomial was used for identifying the parameters by employing aprocedure to be explained below. Naturally, as is the case with any such approximate method, there are also limitationsto be discussed in our numerical examples that follow.Our first and benchmark example concerns a synthetic time series, inspired by parameter values that are at theproper ballpark for describing the real-life data of Portugal’s deaths [54]: That is to say, we prescribed a parameter setof all 9 parameters, resulting in the lumped (per our analysis above) coefficients of Table 1. Upon assigning realisticinitial conditions, we ran the differential equations (1)–(6) forward, and produced the time series of the top left panelof Fig. 2. Notice that in this case all data ( 𝐼 ( 𝑡 ) , 𝐴 ( 𝑡 ) , 𝑆 ( 𝑡 ) , 𝐻 ( 𝑡 ) , 𝑅 ( 𝑡 ) , 𝐷 ( 𝑡 )) are available (in this procedure, the onlyerror involved stems from the approximate nature of our numerical integration); from these quantities we constructed 𝐷 (1) , 𝐷 (2) , 𝐷 (3) and 𝐷 (4) , again within the above (small for our high order numerical scheme) error. The first check to A.S. Fokas, J. Cuevas-Maraver and P.G. Kevrekidis:
Preprint submitted to Elsevier
Page 7 of 14 quantitative framework for exploring exit strategies from the COVID-19 lockdown which we subjected our method is that of identifying 𝛽 : if the above time series of 𝐷 ( 𝑡 ) and its derivatives are accurate,using Eq. (7) to obtain 𝛽 ( 𝑡 ) , the corresponding time series should be found to be constant. We observed this was thecase up to the 5th significant digit, as shown in the top middle panel of Fig. 2, in line with the high accuracy of ournumerical integrator.Then, we turned to the inverse problem. Suppose we are only given at first the time series of 𝐷 ( 𝑡 ) and its derivatives.Can we use the time series of 𝑓 ( 𝑡 ) and in particular (practically) the minimization of || 𝑓 || to uniquely identify thereduced, uniquely identifiable parameters within the set from which we started? To provide the answer for this syntheticexample, we initially substituted all of these time series to the left-hand side of Eq. (7) and examined whether theminimization of the time-series stemming from the ODE provided the expected local minimum (which was a prioriknown as a result of the knowledge of all model parameters for this time series). Indeed, we found this to be the case,as per the data presented in Table 1, under “exact” (indicating that this was generated from the exact time series –up tothe approximations of numerical integration– for 𝐷 ( 𝑡 ) and also 𝐷 (1) , 𝐷 (2) , 𝐷 (3) and 𝐷 (4) ). We then explored a grid ofvalues around this central minimum to decide whether it is a local or global minimum. In this grid, we included, for each of the 8 parameters entering the minimization, values which corresponded to the minimum (the minimum willbe denoted as “1” i.e., 1 × the minimum), denoting of the minimum, , and times the minimum; in thisway, we created a tractable “grid” of values within a fairly large range of variation from the ballpark values of interest,and posed the following question: was the expected from our analysis global minimum the one associated with thecoefficients that we started with (the set of 1’s in the language of the multiples considered) or not? The answer withinthis table of = 390625 entries for the present example was found to be in the affirmative.In order to obtain a more refined sense of whether we really obtained the global minimum, we then went on toexplore a finer grid of parameter values. To that effect, we decided to explore (and illustrate) the method on the basisof a considerably more refined Table but of fewer parameters. Through our various explorations of the method, weidentified as the parameters that were most robust under fitting the ones associated with direct rates, namely 𝑅 , 𝑅 , 𝑟 and 𝐹 = 𝑎 + 𝑠 . Also, 𝛼 barely varied due to its being dominated by the contribution of 𝑇 within its definition;similarly with 𝛽 . In that light, we decided to hereafter focus on the role of 𝐶 and 𝐶 that are also, in a sense, some ofthe most crucial parameters of the model in controlling additional infections, the spreading of the epidemic (also perEq. (11)) and, hence, ultimately the deaths that result from the model. They are also the ones that would be cruciallyaffected under the easing of lockdown conditions; hence, in our view these parameters are the most significant ones toexamine in terms of the robustness of the result across such variations. Among the two ( 𝐶 and 𝐶 ), the crucial roleof asymptomatics (due to their mobility within the population) renders 𝐶 as the typically dominant and significantlylarger value. In that light, we restricted our considerations to the interval 𝐶 ∕10 ≤ 𝐶 ≤ 𝐶 ∕3 (unless indicatedotherwise) in what follows below. This was done to avoid combinations of ( 𝐶 , 𝐶 ) that may lead to âĂŸcompetingminimaâĂŹ, where 𝐶 is comparable or larger than 𝐶 ; this would mean that asymptomatics would have a similar orsmaller number of contacts as the (expected to be self-isolating) sick, a feature not in line with the epidemiologicalpremise of the model. Given the above considerations, we decided to consider a refined variation of 𝐶 over a gridof (1∕5 , . , . , … , . , , . , … , . , . , the original value, and similarly of 𝐶 (within the limits of theabove inequality condition) of the above reported minimum of Table 1. The results of this fine sampling variation areshown in the bottom left panel of Fig. 2, as well as in Table 2. These results confirmed that the original identifiedminimum is the global one with the optimal fitness within this parametric grid. Table 2 also compares this fitness bymeans of measuring || 𝑓 || (which is what is also represented in Fig. 2) not only for this global minimum (our (1,1) stategiven by the red spot in this figure and also in the panels that will follow), but also for the âĂŸnext best thingsâĂŹ,which are given by black spots in the figure. The edge of the original minimum as being the global one by some orders of magnitude is clearly evident within this refined grid. This confirmed the benchmark of our expectations forthis essentially exact time series. Indeed, in this case even an extensive search over arbitrary 𝐶 ’s (not conforming tothe above inequality) would still not improve the global minimum result.The next line of attack was to seek a similar optimization approach in this well-curated and smooth time-seriesexample (again the one of the top left of Fig. 2), but now assuming that we are not given both 𝐷 ( 𝑡 ) and all of itsderivatives, but rather only 𝐷 ( 𝑡 ) . Our first approximation of 𝐷 (1) , 𝐷 (2) , 𝐷 (3) and 𝐷 (4) in this smooth case was basedon the numerical method of finite differences. We used this well-known tool to approximate all derivatives. Then,again, we first posed the following question: equipped with these time series, could we find a local minimum of Eq. (7)through minimizing || 𝑓 || in the vicinity of the exact result? The answer was again in the affirmative and is accountedfor in the corresponding “Finite Differences” column of Table 1. It can be seen that despite the approximate nature ofthe method, we still effectively captured all relevant (lumped, reduced) coefficients. Subsequently, we again followed A.S. Fokas, J. Cuevas-Maraver and P.G. Kevrekidis:
Preprint submitted to Elsevier
Page 8 of 14 quantitative framework for exploring exit strategies from the COVID-19 lockdown the sampling procedure developed earlier with the fine grid over the most sensitive and, arguably, most important(towards the easing of the lockdown measures) parameters, namely 𝐶 and 𝐶 . Once again, we found, as seen in themiddle bottom panel of Fig. 2, that the original seed, namely the minimum that was expected, is the one that beatsthe competition in terms of its fitness by being the lowest one. However, the caveats of the approximate nature of themethod (when only 𝐷 is prescribed) are starting to come into light. In particular, it can be seen that next to the red spot,corresponding to the (1 , parametric set, the one that was expected, there are others that are less clearly discerniblethan the previous example in terms of their fitness. In other words, the edge of our expected optimum was significantlydecreased by the approximate nature of the time series: from an edge of orders of magnitude, we went to one of only , as is shown in first row of Table 3. As a way to read the relevant tables (also for the sake of future examples), wenote that the fitness of the global minimum in terms of its || 𝑓 || is . −7 , while the next best thing concerns a 𝐶 that is . times larger and a 𝐶 that is . smaller and yields a || 𝑓 || = 3 .
65 × 10 −5 .This issue comes into sharper focus when we proceeded to the further approximation of the logistic formula (sig-moid). In the case of a realistic time series, the data becomes quite noisy and hence the notion of finite differencesintroduces a considerable amount of this noise into the minimization of || 𝑓 || . Hence, the logistic approach is deemedto be a worthwhile alternative, to consider because of the smoothness of the curve and its derivatives, and the effectivetransformation of the left-hand side of 𝑓 into a well-defined function of 𝑥 (possibly even a polynomial upon a suitableexpansion of the logarithmic term). However, these advantages come with some caveats too. As we can see in thesynthetic time series at hand, we can again perform a minimization with the logistic approximation in the vicinity ofthe original exact parametric set of this example. We indeed found, in that case too, a local minimum reflected underthe sigmoid column of Table 1, which is in good agreement with the previous methods; notice, however, the slightdeviations in the integration constants due to the approximate nature of the logistic model. In this progressively morerealistic case (since, typically, neither a smooth time series for 𝐷 ( 𝑡 ) and all its derivatives, nor even a smooth one for 𝐷 alone will be practically available), we see that the fitness of our global minimum only edges that of the next closestminimum by a smaller amount than in the finite difference case. However, it is crucial to emphasize that in this exam-ple, as well as in the examples to follow, the closest second minimum, 𝐶 that plays the dominant role in determiningthe fate of dynamics of 𝐷 ( 𝑡 ) , is very close to (only . times) that of the Table 1. Hence, even in a foreseeable casewhere such a secondary minimum may eclipse our expected global one, the dynamics of the two cases will be veryclose to each other.As may be anticipated, similar issues arose when we attempted to utilize the logistic approach to real data for the3 examples outlined in the main text: the countries of Portugal (left panels of Fig. 3), Greece (middle panels of Fig. 3)and the autonomous community of Andalusia in Spain (right panels of Fig. 3). In each case, we attempted to fit aportion of the data, the one past the inflection point in order to capture as best as possible the asymptotic number ofdeaths 𝐷 𝐹 which is a crucial aspect of the sigmoid (together with the growth coefficient 𝑘 ). In our experience, 𝑘 isusually the most robust feature of the associated sigmoids; 𝐷 𝐹 is the trickiest one to obtain, especially with time seriessuch as the ones of Fig. 3 which have not clearly asymptoted to a concrete value. Typically, we have found the sigmoidto under-predict 𝐷 𝐹 . We have identified this as one of the significant deficiencies of the sigmoid which, in turn, callsfor improvement of the method, e.g., along the lines of [52, 53]. The associated fits to sigmoids for the 3 examples areshown in the top panels of Fig. 3, while the associated sigmoid parameters are shown in Table 3. Although in theseexamples the use of the sigmoid was sufficient, developing the relevant extension along the lines of e.g. the rationalapproximation [52, 53] may be a natural vein for future studies and may indeed enhance the robustness of the methodproposed herein (especially in the presence of noisier or more sparse time series).The bottom panels of the figure illustrate further what may happen in the case of the sigmoid. Notice that in thecase of these real data for Portugal (left), Greece (middle) and Andalusia (right), we have no a priori knowledge ofwhat the optimum fit parameters are. However, as a comparison measure, we utilized the parameter values found by aconstrained minimization along the lines discussed in [48] (see also [55]). Admittedly, this minimization using Mat-lab’s routine fmincon (for minimizing the distance between true and observed time series stemming from integratingthe differential equations) does not offer any guarantee of a global minimum and indeed may yield a local minimum.Nevertheless, it provides a useful benchmark for the performance of the minimization of || 𝑓 || within Eq. (7) by usingthe logistic approximation. Subsequently, based upon the logistic fit (most notably its 𝑘 and 𝐷 𝐹 ), we attempted toidentify a local minimum through the logistic method in the vicinity of the one obtained by constrained minimization.Indeed, such a minimum was found in each of the 3 cases; and its parametric values are illustrated in Table 4. Then, weattempted again a refined grid of ( 𝐶 , 𝐶 ) conforming to the same constraint as before, except for the case of Greece:there, given the result of the constrained minimization for 𝐶 ≈ 𝐶 ∕3 , we extended our search to a wider parametric A.S. Fokas, J. Cuevas-Maraver and P.G. Kevrekidis:
Preprint submitted to Elsevier
Page 9 of 14 quantitative framework for exploring exit strategies from the COVID-19 lockdown interval of . 𝐶 ≤ 𝐶 ≤ 𝐶 for completeness. Indeed, we checked that our results for this case did not change evenfor a wider parametric search.In all the explored cases, we found the local minimum (the one obtained both by constrained minimization andindependently by the logistic approach and presented in Table 4) to be the global minimum within the interval ofour search. Nevertheless, once again, it is clear that other combinations of 𝐶 and 𝐶 may yield fitnesses close tothat minimum, since the differences under the logistic approximation were always less than one order of magnitude.Naturally, this raises concerns regarding the identification of the global minimum. However, we note the following:(i) In this case, we have checked (something that is possible more generally) the result in the interval of interestunder two separate methods. (ii) More importantly, as noted earlier, the secondary minima identified involve similarresults as regards 𝐶 and hence are expected to have similar dynamical implications. (iii) More refined variants of themethod (e.g., using the rational or birational models [52, 53]) should be able to address the needs of more complex/lessinformative cases than the ones discussed here. In that vein, we believe that we have exposed both the advantages andthe caveats of the proof-of-principle logistic approach and how it can be improved, as needed. Figure 2:
The top left panel shows the synthetic time series of 𝐷 ( 𝑡 ) . The top middle panel shows the numerical computationof 𝛽 illustrating that indeed as expected this first integral of the motion remains invariant over time, up to numericalaccuracy. In the bottom panel the time series of the left-hand side of Eq. (7) is evaluated and its 𝑙 norm is calculated.This is done for the grid of finely sampled values of ( 𝐶 , 𝐶 ) . The red dot indicates the local optimum associated withthe known parameter values from which the time series was initially constructed. In the left panel the rest of the dotscorrespond to the ones reported in Table 2 (the next best fits). The bottom panels correspond to the fitness indicator || 𝑓 || for the exact time series (including its derivatives) on the left, for the time series and then its derivatives from finitedifferences (middle), and from a logistic approximation (right). Notice how the identified local minimum of the red dotturns out to also be the global minimum within the sampling interval considered.
4. Conclusions and Future Recommendations
We conclude with several remarks, summarizing our work and offering a number of future recommendations andinteresting research directions:1. We have introduced a rigorous and computationally tractable methodology for computing the impact on thenumber of deaths of increasing the parameter 𝑐 in our SIR-type model. This parameter is proportional tothe number of contacts between asymptomatic individuals infected with SARS-CoV-2 and those susceptible.Remarkably, the only data needed for our algorithm to be implemented is the cumulative number of deaths duringthe lockdown period. Our considerations have been based on an extended version of an SIR model (incorporating A.S. Fokas, J. Cuevas-Maraver and P.G. Kevrekidis:
Preprint submitted to Elsevier
Page 10 of 14 quantitative framework for exploring exit strategies from the COVID-19 lockdownSeries 𝐶 𝐶 || 𝑓 || || 𝑓 || (global)Synthetic (Finite differences) . . .
65 × 10 −5 .
30 × 10 −7 Synthetic (Sigmoid) . . .
89 × 10 −7 .
91 × 10 −9 Real series (Portugal) . . .
30 × 10 −6 .
03 × 10 −6 Real series (Greece) . . .
74 × 10 −11 .
21 × 10 −11
Real series (Andalusia) . . .
28 × 10 −7 .
91 × 10 −8 Table 3
Similar to Tab. 2 but now for each of the remaining 5 examples (corresponding to the 5 rows). Here the fitness of theglobal optimum is shown in the 5th and last column, while the next best thing (in multiples of the optimal ( 𝐶 , 𝐶 ) ) andits associated fitness are shown in columns (2 , and , respectively.Parameter Portugal Greece Andalusia 𝐷 𝑓 . .
93 1404 . 𝑘 . . . 𝑃 . . . 𝐶 . −5 . −5 . −6 𝐶 . −6 . −6 . −6 𝑅 . . . 𝑅 . . . 𝐹 . . . 𝛼 . . . 𝑟 . . . 𝛽 −1 . . . Table 4
Similarly to Table 1 (for Fig. 2), but now for the 3 sigmoid-based country/autonomous province examples of Fig. 3, thelumped parameter combinations of the model are given, together with the coefficients of optimal sigmoid fit performedwithin the relevant figure. asymptomatically infected, hospitalized individuals, as well as deaths). Using an inverse problem perspective,we have formulated a single ODE for the fatalities, which we consider as the “ground truth” within the availabledata. Via an optimization scheme, we have utilized this 𝐷 ( 𝑡 ) (and its derivatives, obtained, e.g., via a logisticfitting to the data) to obtain the uniquely identifiable parametric combinations within the model.2. In addition to providing synthetic examples to illustrate the method, we have also utilized real data in orderto showcase the practical implementation of the proposed approach. The application of our algorithm towardsthe study of the pandemic spread to Greece, Portugal and Andalusia shows that upon easing of the lockdownmeasures by increasing the value of 𝑐 (and hence effectively the number of contacts) by a multiple larger than2, the impact on the number of deaths will be devastating. However, an increase by 𝜁 = 2 has a more modesteffect, especially in the case of Greece which is found to have a low transmission rate.3. Importantly, a similar analysis can be applied to sub-populations. Let us consider two such subpopulations: oneconsisting of individuals below the age of 40, that will be referred to as young, and one consisting of individualsabove 40, that will be denoted as older. If there exist reliable deaths data for these two subpopulations, then itis possible to compute separately the effect of changing the interactions among the young, as well as among theolder with the older and the older with the young. We did, for example, look at these two subpopulations in thecase of Greece [48]. In this case, we supplemented the data of deaths with data on the cumulative number ofreported infected, because the data on deaths, especially for the young, were sparse. These results provided hope(and a potential future recommendation, namely) that it may be possible to accelerate the easing of the lockdownfor the young: the increase of interactions among the young (as opposed to the interactions involving the older)led to only a modest increase in the number of deaths. However, it should be noted that, as discussed above,our conclusions for single populations are based on rigorous mathematical results and reliable data; since thisis less transparently so the case for the data (e.g., for the reported infected) regarding the two sub-populationsof Greece, the latter results can only be considered as suggestive. From a mathematical perspective, completing A.S. Fokas, J. Cuevas-Maraver and P.G. Kevrekidis:
Preprint submitted to Elsevier
Page 11 of 14 quantitative framework for exploring exit strategies from the COVID-19 lockdown
Figure 3:
The top 3 panels present the data for 𝐷 ( 𝑡 ) for the countries of Portugal (left column), Greece (middle column),and the autonomous community of Andalusia in Spain (right column). Also shown is a logistic fit of the central portionof the associated data aimed at capturing the logistic coefficients 𝑘 and 𝐷 𝐹 . Then in the bottom panel, similarly to thebottom panels of Fig. 2, but now only for the logistic approximation of the top panel, the optimum within a refined gridof values of the 𝐶 and 𝐶 parameters is traced; the corresponding seed at the center of our consideration is representedby the parametric set found by constrained minimization and shown by a red dot in each case. and practically applying the rigorous analysis in the same spirit as above for two or more populations is certainlya challenging problem for future studies.4. Our algorithm can be used in any country with reliable death data during the lockdown period; if death dataalso exist for sub-populations, the multicomponent version of our algorithm should also be applicable. It wouldconstitute an interesting future recommendation to apply our model widely to a number of other case examples,in order to identify the parametric differences between countries with high and ones with lower fatality casenumbers.5. A central mathematical finding of our analysis concerns the rigorous illustration of the fact that the solution ofthe inverse problem (based on death data) does not yield in a unique way all the model parameters, thus it isnot directly possible to compute the effect of the easing of lockdown on the number of hospitalizations and onthe total number of infected individuals. However, we can compute the arguably most important feature of theepidemic, upon release of lockdown, namely the number of anticipated future deaths on the basis of the model.6. At the moment, we cannot provide a specific relation between concrete protective measures (such as hygieneconditions, partial social distancing, and perhaps especially wearing a mask) and the value of 𝜁 . Such a relationexists, as 𝑐 does not solely reflect the number of contacts but also the probability of infection given a contactwhich is proportional to the viral load (i.e., the viral concentration in the respiratory-tract fluid) of expelled res-piratory droplets [50]. Clearly, the above protective measures reduce the viral load and therefore 𝑐 . Hopefully,such a relation can be established by designing specific experiments in the near future. This is a future rec-ommendation that would be especially valuable towards assessing the impact of models within the broad classconsidered herein.7. Our analysis makes the crucial assumption that the basic characteristics of the virus in the post-lockdown stageare identical with those during the lockdown. If, for example, a mutation takes place which makes the virusweaker, then the effect of easing the lockdown measures will not be as dramatic as predicted in our work. Acknolwegdments.
ASF acknowledges support from EPSRC in the form of a a Senior Fellowship. PGK andJCM greatly appreciate numerous helpful discussions with Y. Drossinos, Z. Rapti, and G.A. Kevrekidis. This material
A.S. Fokas, J. Cuevas-Maraver and P.G. Kevrekidis:
Preprint submitted to Elsevier
Page 12 of 14 quantitative framework for exploring exit strategies from the COVID-19 lockdown is based upon work supported by the US National Science Foundation under Grants No. PHY-1602994 and DMS-1809074 (PGK). PGK also acknowledges support from the Leverhulme Trust via a Visiting Fellowship and thanks theMathematical Institute of the University of Oxford for its hospitality during part of this work. The work presented hereis part of a project that we begun jointly with Nikos Dikaios, with whom we hope to present elsewhere other aspectsof this project.
References [1] Z. Sedaghat, N. Karimi, Guillain Barre syndrome associated with COVID-19 infection: A case report, J. Clin. Neurosci.10.1016/j.jocn.2020.04.062 (2020).[2] M.C. Dalakas, COVID19 and GuillainâĂŞBarré Syndromes: the first documented COVID19-triggered autoimmune Neurological disease. toappear in Neurology.[3] V.J. Munster, M. Koopmans, N. van Doremalen, D. van Riel, E. de Wit, A Novel Coronavirus Emerging in China – Key Questions for ImpactAssessment, N. Eng. J. Med. , 692 (2020).[4] See, e.g., https://elifesciences.org/articles/57309 [5] See, e.g., .[6] N. Zhu, D. Zhang, W. Wang et al. , A Novel Coronavirus from Patients with Pneumonia in China, N. Engl. J. Med. , 727 (2020).[7] M. Vaduganathan, O. Vardeny, T. Michel, et al. , ReninâĂŞAngiotensinâĂŞAldosterone System Inhibitors in Patients with Covid-19. N. Engl.J. Med. , 1653 (2020).[8] J. Fantini, C. Di Scala, H. Chahinian, N. Yahi. Structural and molecular modeling studies reveal a new mechanism of action of chloroquineand hydroxychloroquine against SARS-CoV-2 infection, Int. J. Antimicrob. Agents (2020). https://doi.org/10.1016/j.ijantimicag.2020.105960 [9] B. Cao et al. , A Trial of LopinavirâĂŞRitonavir in Adults Hospitalized with Severe Covid-19, N. Engl. J. Med. , 1787 (2020). https://doi.org/10.1056/NEJMoa2001282 [10] J. Magagnoli et al. , Outcomes of hydroxychloroquine usage in United States veterans hospitalized with Covid-19, medRxiv https://doi.org/10.1101/2020.04.16.20065920 [11] M. Hoffmann et al. , Chloroquine does not inhibit infection of human lung cells with SARS-CoV-2. Nature (2020). https://doi.org/10.1038/s41586-020-2575-3 [12] P. Maisonnasse et al. , Hydroxychloroquine use against SARS-CoV-2 infection in non-human primates. Nature (2020). https://doi.org/10.1038/s41586-020-2558-4 [13] T. Chivese et al. , A meta-review of systematic reviews and an updated meta-analysis on the efficacy of chloroquine and hydroxychloroquinein treating COVID19 infection. https://doi.org/10.1101/2020.07.28.20164012 [14] N.M. Ferguson et al. , Impact of non-pharmaceutical interventions (NPIs) to reduceCOVID-19 mortality and healthcare demand. ImperialCollege, London, Mar. 2020. [15] S. Flaxman et al. , arXiv:2004.11342.[16] https://github.com/CSSEGISandData/COVID-19/tree/master/csse_covid_19_data[17] K.K. Ali, C. Cattani, J.F. Gómez-Aguilar, D. Baleanu, M.S. Osman, Analytical and numerical study of the DNA dynamics arising in oscillator-chain of Peyrard-Bishop model, Chaos, Solitons & Fractals , 110089 (2020).[18] B. Inan, M.S. Osman, T. Ak, D. Baleanu, Analytical and numerical solutions of mathematical biology models: The NewellâĂŘWhitehead-âĂŘSegel and AllenâĂŘCahn equations, Math. Meth. Appl. Sci. , 2588 (2020).[19] D. Lu, M.S. Osman, M.M.A. Khater, R.A.M. Attia, D. Baleanu, Analytical and numerical simulations for the kinetics of phase separation iniron (FeâĂŞCrâĂŞX (X=Mo, Cu)) based on ternary alloys, Physica A , 122634 (2020).[20] H. Rezazadeh, M.S. Osman, M. Eslami, M. Mirzazadeh, Q. Zhou, S.A. Badri and A. Korkmaz, Hyperbolic rational solutions to a variety ofconformable fractional Boussinesq-Like equations, Nonlin. Engin. , 224 (2018).[21] M.S. Osman, H. Rezazadeh and M. Eslami, Traveling wave solutions for (3+1) dimensional conformable fractional Zakharov-Kuznetsovequation with power law nonlinearity, Nonlin. Engin. , 559 (2019).[22] O.A. Arqub, M.S. Osman, A.-H. Abdel-Aty, A.-B. A. Mohamed, S. Momani, A Numerical Algorithm for the Solutions of ABC SingularLaneâĂŞEmden Type Models Arising in Astrophysics Using Reproducing Kernel Discretization Method, Mathematics , 923 (2020).[23] E. Kaxiras, G. Neofotistos, E. Angelaki, The first 100 days: modeling the evolution of the COVID-19 pandemic, Chaos, Solitons & Fractals,in press (2020).[24] E. Kaxiras, G. Neofotistos, Multiple epidemic wave model of the Covid-19 pandemic, J. Med. Internet Res.[25] E. Kaxiras, G. Neofotistos, Modeling the COVID-19 Pandemic Response of the US States, (preprint).[26] K. Mizumoto, K. Kayaga, A. Zarebski, and G. Chowell, Estimating the Asymptomatic Proportion of 2019 Novel Coronavirus onboard thePrincess Cruises Ship, 2020, Eurosurveillance doi: 10.2807/1560-7917.ES.2020.25.10.2000180[27] P. Yang, J. Qi, S. Zhang, X. Wang, G. Bi, Y. Yang, and B. Sheng, Feasibility Study of Mitigation and Suppression Intervention Strategies forControlling COVID-19 Outbreaks in London and Wuhan, https://doi.org/10.1101/2020.04.01.20043794 [28] L. Danon, E. Brooks-Pollock, M. Bailey, and M. Keeling, A spatial model of CoVID-19 transmission in England and Wales: early spread andpeak timing, [29] Ka-Ming Tam, Nicholas Walker,and Juana Moreno, Projected Development of COVID-19 in Louisiana, arXiv:2004.02859.[30] A. Arenas et al. , A mathematical model for the spatiotemporal epidemic spreading of COVID19, https://doi.org/10.1101/2020.03.21.20040022 [31] C. Qi, D. Karlsson, K. Sallmen, and R. Wyss, Model studies on the COVID-19 pandemic in Sweden, arXiv:2004.01575 A.S. Fokas, J. Cuevas-Maraver and P.G. Kevrekidis:
Preprint submitted to Elsevier
Page 13 of 14 quantitative framework for exploring exit strategies from the COVID-19 lockdown [32] G.D. Barmparis and G.P. Tsironis, Estimating the infection horizon of COVID-19 in eight countries with a data-driven approach, Chaos,Solitons & Fractals , 109842 (2020).[33] E. Gjini, Modeling Covid-19 dynamics for real-time estimates and projections: an application to Albanian data, https://doi.org/10.1101/2020.03.20.20038141 [34] S.B. Bastos and D.O. Caljueiro, Modeling and forecasting the early evolution of the Covid-19 pandemic in Brazil, arXiv:2003.14288.[35] V.Z.Marmarelis, Predictive modelling of Covid-19 data in the US: Adaptive phase-space approach, IEEE Open J. Eng. Med. Biol. (to appear).[36] C. Yang and J. Wang, A mathematical model for the novel coronavirus epidemic in Wuhan, China, Math. Biosci. Eng. , 2708 (2020).[37] L. Russo, C. Anastassopoulou, A. Tsakris, G.N. Bifulco, E.F. Campana, G. Toraldo, and C. Siettos, Tracing DAY-ZERO and Forecastingthe COVID-19 Outbreak in Lombardy, Italy: A Compartmental Modelling and Numerical Optimization Approach, https://doi.org/10.1101/2020.03.17.20037689 ; C. Anastassopoulou, L. Russo, A. Tsakris, C. Siettos, Data-Based Analysis, Modelling and Forecasting ofthe COVID-19 outbreak, PLoS ONE , e0230405 (2020).[38] W. Guan, et al. , Clinical Characteristics of Coronavirus Disease 2019 in China, NEJM , 1708 (2020).[39] T. M. McMichael et al. , Epidemiology of Covid-19 in a Long-Term Care Facility in King County, Washington, N. Engl. J. Med. , 2005,(2020).[40] N Lurie, M. Saville, R. Hatchett, J. Halton, Developing Covid-19 Vaccines at Pandemic Speed, N. Engl. J. Med. , 1969 (2020).[41] https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology [42] M. M. Arons et al. , 84 (2013).[45] A.S. Fokas, Electro-Magneto-Encephalography for the Three-Shell Model: Distributed Current in Arbitrary, Spherical and Ellipsoidal Ge-ometries, J. R. Soc. Interface , 479 (2009).[46] A. Holiday, M. Kooshkbaghi, J.M. Bello-Rivas, C.W. Gear, A. Zagaris, I.G. Kevrekidis, Manifold learning for parameter reduction, J. Comp.Phys. , 419 (2019).[47] See, e.g.,: https://time.com/5824836/greece-coronavirus/ and also [48] A.S. Fokas, J. Cuevas-Maraver, P.G. Kevrekidis, Two alternative scenarios for easing COVID-19 lockdown measures: one reasonable and onecatastrophic, https://doi.org/10.1101/2020.05.08.20095380.[49] T. Armangue, M. Spatola, A. Vlagea, et al. , Frequency, symptoms, risk factors, and outcomes of autoimmune encephalitis after herpes simplexencephalitis: a prospective observational study and retrospective analysis. Lancet Neurol. , 760 (2018).[50] Y. Drossinos and N.I. Stilianakis, Aerosol Sci. Technol. , 639 (2020), https://doi.org/10.1080/02786826.2020.1751055 .[51] F.R. Gantmacher, Applications of the Theory of Matrices , Wiley (New York, 1959).[52] A.S. Fokas, N. Dikaios, G.A. Kastis, Mathematical models and Deep Learning for predicting the number of individuals reported to be infectedwith SARS-CoV-2,[53] A.S. Fokas, N. Dikaios, G.A. Kastis. COVID-19: Predictive Mathematical Models for the Number of Deaths (preprint).[54] https://github.com/dssg-pt/covid19pt-data[55] P.G. Kevrekidis, J. Cuevas-Maraver, Y. Drossinos, Z. Rapti, G.A. Kevrekidis, Spatial Modeling of COVID-19: Greece and Andalusia as CaseExamples, arXiv:2005.04527.
A.S. Fokas, J. Cuevas-Maraver and P.G. Kevrekidis: