Abstract
We formulated a mathematical model considering young (below 60 years old) and elder (above 60 years) subpopulations to describe the introduction and dissemination of new coronavirus epidemics in the São Paulo State, Brazil. From the data collected in São Paulo State, we estimated the model parameters and calculated the basic reproduction number as R0 = 6.828. Considering isolation as a control mechanism, we varied the releasing proportions of young and elder persons to assess their epidemiological impacts. The best scenarios were release of young persons, but maintaining elder persons isolated. To avoid the collapse of the health care system, the isolation must be at least 80%.
1 Introduction
Mathematical models allow us to understand the progression of viral infections if the natural history of the disease is well documented. This understanding, in turn, permits forecasting epidemiological scenarios when control mechanisms are introduced aiming to reduce or eliminate the infection. In the case of coronavirus disease 2019 (CoViD-19), which is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), a strain of the RNA-based SARS-CoV-1, the possibility of obtaining scenarios is of fundamental importance. The reason is that in serious cases due to SARS-CoV-2 (new coronavirus) infection, immune cells overreact and attack the lung cell causing acute respiratory disease syndrome and possibly death. In general, the fatality rate in elder patients (60 years or more) is much higher than the average, and under 40 years seems to be around 0.2%. [10].
Due to the rapid spreading out of new coronavirus, its distribution is currently worldwide (pandemic). This virus can be transmitted by droplets that escape the lungs through coughing or sneezing and infect humans (direct transmission), or they are deposited in surfaces and infect humans when in contact with this contaminated surface (indirect transmission). The virus enters into susceptible persons through the nose, mouth, or eyes, and infects cells in the respiratory tract, being capable to release millions of new viruses. Like all RNA-based viruses, the new coronavirus tends to mutate faster than DNA-viruses, but lower than influenza viruses
Currently, there is not a vaccine, neither efficient treatment, although many drugs (chloroquine, for instance) are under clinical trial. Hence, isolation is the main, if not unique, way of controlling the dissemination of this virus in a population aiming the change in the natural history of disease propagation (this change is commonly known as the flattening curve of epidemics). Nevertheless, this isolation arises an important question: are there reliable strategies to release these isolated persons aiming to avoid the retaken of its original progression of infection?
Many mathematical and computational models are being used to describe current new coronavirus pandemics. In the mathematical models, there is a fundamental threshold (see [2]) called the basic reproduction number, which is defined as the secondary cases produced by one case introduced in a completely susceptible population, which is denoted by R0. When a control mechanism is introduced, this number is reduced and is called as the reduced reproduction number Rr. Ferguson et al. [5] proposed a model to investigate the effects of isolation of susceptible persons. They analyzed two scenarios, called by them as mitigation and suppression, and predicted the numbers of severe cases and deaths due to CoViD-19 without control mechanism, and compared them with those numbers when control is introduced.
In this paper, we formulate a mathematical model based on ordinary differential equations, aiming the understanding of the dynamics (or trajectories of dissemination) of new coronavirus transmission, and applying this model to forecast changes in the dynamics under intervention. The model considers pulse isolation and a series of pulses of release (see [11] for a series of pulses vaccination). The model aims to describe the onset and subsequent spread of the new coronavirus in São Paulo State, Brazil. The understanding of the dynamics means that the model parameters are fitted against data collected in São Paulo State, Brazil. These estimated parameters are then used to study the potential scenarios when an intervention is introduced as a control mechanism. São Paulo State adopted the isolation of persons as the controlling mechanism, which can not be isolated indefinitely. Our main aim is to obtain epidemiological scenarios when releasing strategies will be implemented after isolation.
The paper is structured as follows. In Section 2, we introduce a model, which is numerically studied in Section 3. Discussions are presented in Section 4, and conclusions in Section 5.
2 Material and methods
In a community where SARS-CoV-2 (new coronavirus) is circulating, the risk of infection is greater in elder than young persons, as well as under increased probability of being symptomatic and higher CoViD-19 induced mortality. Hence, the community is divided into two groups, composed by young (under 60 years old, denoted by subscript y), and elder (above 60 years old, denoted by subscript o) persons. The vital dynamics of this community are described by per-capita rates of birth (ϕ) and mortality (µ).
For each sub-population j (j = y, o), all persons are divided into seven classes: susceptible Sj, susceptible persons who are isolated Qj, exposed Ej, asymptomatic Aj, asymptomatic persons who are caught by test and then isolated Q1j, symptomatic persons at initial phase of CoViD-19 (or pre-diseased) D1j, pre-diseased persons caught by test and then isolated, plus mild CoViD-19 (or non-hospitalized) Q2j, and symptomatic persons with severe CoViD-19 (hospitalized) D2j. However, all young and elder persons in classes Aj, Q2j, and D2j enter into the same immune class I.
The natural history of new coronavirus infection is the same for young (j = y) and elder (j = o) subpopulations. We assume that only persons in the asymptomatic (Aj) and pre-diseased (D1j) classes are transmitting the virus, and other infected classes (Q1j, Q2j and D2j) are under voluntary or forced isolation. Susceptible persons are infected according to λjSj (known as mass action law [2]) and enter into classes Ej, where λj is the per-capita incidence rate (or force of infection) defined by λj = λ (δjy + ψδjo), with λ being
where δij is Kronecker delta, with δij = 1 if i = j, and 0, if i ≠ j; and β1j and β2j are the transmission rates, that is, the rates at which a virus encounters a susceptible people and infects. After an average period 1/σj in class Ej, where σj is the incubation rate, exposed persons enter into the asymptomatic Aj (with probability pj) or pre-diseased D1j (with probability 1 − pj) classes. After an average period 1/γj in class Aj, where γj is the infection rate of asymptomatic persons, symptomatic persons acquire immunity (recovered) and enter into immune class I. Another route of exit from class Aj is being caught by a test at a rate ηj and enters into class Q1j and, then, after a period 1/γj, enters into class I. Possibly asymptomatic persons are in voluntary isolation, which is described by the voluntary isolation rate χj. With respect to symptomatic persons, after an average period 1/γ1j in class D1j, where γ1j is the infection rate of pre-diseased persons, pre-diseased persons enter into non-hospitalized Q2j (with probability mj) or hospitalized D2j (with probability 1 − mj) class. Hospitalized persons acquire immunity after a period 1/γ2j, where γ2j is the recovery rate of severe CoViD-19, and enter into immune class I or die under the disease induced (additional) mortality rate αj. Another route of exiting D2j is by treatment, described by the treatment rate θj. After an average period 1/γj in class Q2j, non-hospitalized persons acquire immunity and enter into immune class I, or enter into class D2j at a relapsing rate ξj.
In the model, we consider pulse isolation and intermittent (series of pulses) release of persons. We assume that there is a unique pulse in isolation at time 
, described by 
, but there are m intermittent releases described by 
, where 
, and δ (x) is Dirac delta function, that is, δ (x)= ∞, if x = 0, otherwise, δ (x) = 0, with 
. The parameters kj and lij, i = 1, 2, ···,m, are the fractions of i-th release of isolated persons, and τwj is the period between successive releases.
Figure 1 shows the flowchart of the new coronavirus transmission model.
The flowchart of new coronavirus transmission model with variables and parameters.
The new coronavirus transmission model, based on above descriptions summarized in Figure 1, is described by system of ordinary differential equations, with j = y, o. Equations for susceptible persons are
for infectious persons,
and for immune persons,
with the initial number of population at t = 0 being N(0) = N0 = N0y + N0o, where N0y and N0o are the number of young and elder persons at t = 0. If ϕ = µ +(αyD2y + αoD2o) /N, the total size of the population is constant.
Table 1 summarizes the model variables.
The non-autonomous system of equations (2), (3), and (4) is simulated letting intermittent interventions to the initial and boundary conditions. Hence, the equations for susceptible and isolated persons become
j = y, o, and other equations are the same.
For the system of equations (6), (3), and (4), the initial conditions (at t = 0) are, for j = y, o,
and 
is a non-negative number. For instance, 
means that there is not any exposed individual (young and elder) at the beginning of epidemics. We split the boundary conditions into isolation and release, and assume that 
and τi = τiy = τio, for i = 1, 2, ···, m, then 
A unique isolation at t = τis is described by the boundary conditions
plus
where we have 
, and 
. The boundary conditions for a series of pulses released at 
, for i = 1, 2, ···, m, are
plus
If τ = τi, then ti = τis+ iτ. If isolation is applied to a completely susceptible population, at t = 0, there are not any infectious person, so S(0) = N0. If isolation is done at 
without screening of persons harboring the virus, then many of them could be isolated with susceptible persons.
Table 2 summarizes the model parameters and values (for elder classes, values are between parentheses).
Summary of the model parameters (j = y, o) and values (rates in days−1, time in days and proportions are dimensionless). Some values are calculated (&), or varied (#), or assumed (*), or estimated (**) or not available (***).
From the system of equations (2), (3), and (4) we can derive some epidemiological parameters: new cases, new CoViD-19 cases, severe CoViD-19 cases, number of deaths due to CoViD-19, isolated persons, and number of occupied beds in hospitals.
The number of persons infected with new coronavirus is given by Ey + Ay + Q1y + D1y + Q2y + D2y for young persons, and Eo + Ao + Q1o + D1o + Q2o + D2o for elder persons. The incidence rates are
where the per-capita incidence rate λ is given by equation (1), and the numbers of new cases Cy and Co are
with Cy(0) = 0 and Co(0) = 0
The number of CoViD-19 cases Δy and Δo are given by exist from Ay, D1y, Ao, and D1o, that is,
with Δy(0) = 0 and Δo(0) = 0, which are entering in classes Q1y, D2y, Q2y, Q1o, D2o, and Q2o.
The numbers of severe CoViD-19 (hospitalized) cases Ωy and Ωo are given by exits from D1y, Q2o, D2o, and Q2y, that is,
with Ωy(0) = Ωo(0) = 0, which are entering in classes D2y and D2o at each day.
The number of deaths caused by severe CoViD-19 cases Π can be calculated from hospitalized cases. This number of deaths is
with Πy(0) = Πo(0) = 0.
The number of susceptible persons in isolation in the absence of release is obtained from
where the corresponding fractions of isolated susceptible persons are 
and 
Finally, the number of beds B occupied during the evolving of epidemics is, for j = y, o,
for beds in hospitals, and
for beds in ICU. The fraction of severe CoViD-19 needing ICU is hj, and 1/ς1j and 1/ς2j are the average occupying time of beds in hospital and ICU, equal for young and elder persons, where ς1j and ς2j are the discharging rates from hospital and ICU, for j = y, o. The fraction h1j is severe CoViD-19 needing prolonged hospital care. The total number of occupied beds is B = B1 + B2.
The system of equations (2), (3), and (4) is non-autonomous. Nevertheless, the fractions of persons in each compartment approach the steady state (see Appendix A). Hence, at t = 0, the basic reproduction number R0 is obtained substituting 
and 
by 
and 
in equation (A.4), resulting in
using equations (A.8) and (A.9).
Let us use the approximated effective reproduction number Ref given by equation (A.11), that is, Ref = R0 (Sy/N + So/N). For t> 0, we have the effective reproduction number Ref, with Ref (0) = R0 at t = 0, which decreases as susceptible persons decrease. However, at t = τis a pulse in isolation is introduced, hence we have Ref (τis+) = Rr, where the reduced reproduction number Rr is given by
and Sy (τis−) and So (τis−) are the number of susceptible young and elder persons at the time just before the introduction of isolation. Notice that at t = τis, Ref (τis−) jumped down to Ref (τis+). At i-th release time ti, we have 
, with the reduced reproduction number Rr(i) being given by
where 
jumped up to 
at t = ti. After t>tm, there is not release anymore, and Ref = 1 when t →∞, and the new coronavirus returns to the original dynamics driven by R0.
Given N and R0, let us evaluate the number of susceptible persons to trigger and maintain epidemics. Let us assume that all model parameters for young and elder classes and all transmission rates are equal, then R0 = σβ/ [(σ + ϕ)(γ + ϕ)] and Ref = R0S/N, using the approximated Ref given by equation (A.11). Letting Ref = 1, the critical number of susceptible persons Sth at equilibrium is
If S >Sth, epidemics occurs and persists (Ref > 1, non-trivial equilibrium point P*), and the fraction of susceptible individuals is s* = 1/R0, where 
; but if S < Sth, epidemics occurs but fades out (Ref < 1, trivial equilibrium point P0), and the fractions of susceptible individuals sy and so at equilibrium are given by equation (A.4). Remember that the threshold Sth is valid only if the pulse isolation introduced at t = 0 is maintained forever.
We apply the above results to study the introduction and establishment of the new coronavirus in São Paulo State, Brazil.
3 Results
The results obtained in the foregoing section are applied to describe the new coronavirus infections in São Paulo State, Brazil. The first confirmed case of CoViD-19, occurred on February 26, 2020, was a traveler returning from Italy on February 21, and being hospitalized on February 24. The first death due to CoViD-19 was a 62 years old male with comorbidity who never travelled to abroad, hence considered as autochthonous transmission. He manifested the first symptoms on March 10, was hospitalized on March 14, and died on March 16. On March 24, the São Paulo State authorities ordered isolation of persons acting in non-essential activities, as well as students of all levels until April 6, further isolation was extended to April 22, and postponed to May 10.
The dynamics of the new coronavirus propagation is obtained by evaluating the system of equations (2), (3), and (4) numerically using the 4th order Runge-Kutta method. Let us determine the initial conditions supplied to this system. In the São Paulo State, the number of inhabitants is N (0) = N0 = 44.6×106 according to SEADE [8]. The value of parameter φ given in Table 1 was calculated by the equation (A.4), φ = bϕ/ (1 − b), where b is the proportion of elder persons. Using b = 0.153 in the São Paulo State [8], we obtained φ = 6.7 × 10−6 days−1, hence, 
and 
. The initial conditions for susceptible persons are set to be Sy (0) = N0y and So (0) = N0o. For other variables, from Table 2, using py = 0.8 and my = 0.8, the ratio asymptomatic:symptomatic is 4 : 1, and the ratio mild:severe (non-hospitalized:hospitalized) CoViD-19 is 4 : 1. We also use these ratios for elder persons, even po and mo are slightly lower. Hence, if we assume that there is one person in D2j (the first confirmed case), then there are 4 persons in Q2j. The sum (5) is the number of persons in class D1j, implying that there are 20 in class Aj, hence, the sum (25) is the number of persons in class Ej. Finally, we suppose that no one is isolated or tested, and immunized. (Probably the first confirmed COViD-19 person transmitted virus (since February 21 when returned infected from Italy), as well as other asymptomatic travelers returning from abroad.)
Therefore, the initial conditions supplied to the dynamic system (2), (3), and (4) are
where the initial simulation time t = 0 corresponds to the calendar time February 26, 2020, when the first case was confirmed.
This section presents parameters estimation and epidemiological scenarios considering isolation as the control mechanism. In estimation and epidemiological scenarios, we assume that all transmission rates in young persons are equal, as well as in elder persons, that is, we assume that
hence, the forces of infection are λy = (Ay + D1y + Ao + D1o) βy/N and λo = ψλy. The reason to include factor ψ is the reduced capacity of defense mechanism by elder persons (physical barrier, innate and adaptive immune responses, etc.). The force of infection takes into account all virus released by infectious individuals (Ay, D1y, Ao, and D1o), the rate of encounter with susceptible persons, and the capacity to infect them (see [12] [13]). Additionally, the amount aspired by susceptible persons can be determinant in the chance of infection and in the prognostic of CoViD-19 [19].
From data collected in the São Paulo State from February 26 until April 5, 2020, we fit transmission (βy and βo) and additional mortality rates (αy and αo), and the proportions of isolated persons (ky and ko) are estimated using data from March 24 until April 21. Once determined these parameters, we study potential scenarios introducing isolation as control mechanisms.
3.1 Parameters estimation
Reliable estimations of both transmission and additional mortality rates is crucial, aiming the prediction of new cases (to an adequate number of beds in a hospital, for instance) and deaths. When the estimation is based on a few number of data, that is, at the beginning of epidemics, some cautions must be taken, because the rates maybe over or under estimated. The reason is that in the very beginning phase of epidemic, the spreading out of infections and death increase exponentially without bound.
Currently, there is not a sufficient number of kits to detect infection by the new coronavirus. For this reason, test to confirm infection by this virus are done only in hospitalized persons, and also in persons who died, manifesting symptoms of CoViD-19. Hence, we have only data of hospitalized persons (D2y and D2o) and those who died (Πy and Πo). Taking into account hospitalized persons with CoViD-19, we fit the transmission rates, and for persons died due to CoViD-19, we fit the additional mortality rates. These rates are fitted applying the least square method (see [7]), that is
where min stands for the minimum value, n is the number of observations, ti is i-th observation time, Zj stands for Ωj given by equation (13), and Πj is given by equation (14); and 
stands for the observed number of hospitalized persons 
and number of died persons 
The fitted parameters are those minimizing the sum of squared differences.
Instead of using equation (20), the least square estimation method, we vary the parameters and choose better fitting by evaluating the sum of squared distances between curve and data.
3.1.1 Fitting the transmission rates
The introduction of quarantine was at t = 27, corresponding to the calendar time March 24, but the effects are expected to appear later. Hence, we will estimate taking into account the confirmed cases from February 26 (t = 0) to April 5 (t = 39), hence n = 40 observations. It is expected that at around simulation time t = 36 (April 2), the effects of isolation will appear (the sum of incubation and pre-diseased infection periods (see Table 2) is around 9 days).
To estimate the transmission rates βy and βo, we let αy = αo = 0 and the system of equations (2), (3), and (4), with initial conditions given by equation (7), is evaluated, and we calculate
by varying βy and βo. We chose the transmission rates minimizing the sum of differences. Letting additional mortality rates equal to zero (αy = αo = 0), we estimate βy and βo = ψβy, against hospitalized CoViD-19 cases (D2 = D2y + D2o) data from the São Paulo State. The estimated values are βy = 0.75 and βo = 0.8775 (days−1), where Ψ = 1.17, resulting in the basic reproduction number R0 = 6.828 (partials R0y = 5.587 and R0o = 1.241), according to equation (A.8). Figure 2(a) shows the estimated curve of D2 and observed data, plus two curves with lower transmission rates: βy = 0.55 and βo = 0.6435 (days−1), with R0 = 5.008; and βy = 0.45 and βo = 0.5265 (days−1), with R0 = 4.097 (partials R0y = 3.352 and R0o = 0.745). Figure 2(b) shows curves, from t = 30 until 180 of D2 for young, elder, and total persons for R0 = 6.828 and 4.097.
The estimated curve of severe CoViD-19 cases D2 (βy = 0.75 and βo = 0.8775 (days−1), with R0 = 6.828) and observed data, plus two curves with lower transmission rates: βy = 0.55 and βo = 0.6435 (days−1), with R0 = 5.008; and βy = 0.45 and βo = 0.5265 (days−1), with R0 = 4.097. Estimation was done using data from t = 0 to t = 39 (a), from top to bottom as R0 increases, and the extended curves until t = 250 (b), from top to bottom for total, young and elder persons, with continuous (R0 = 6.828) and dashed (R0 = 4.097) curves.
From Figure 2(a), as R0 decreases, the estimations become worst. From Figure 2(b), the peak of higher R0 occurs at t = 68, while for lower R0, delayed in 33 days, at t = 101. The peaks for R0 = 6.828 and 4.097 (between parentheses) are 6.86 × 105 (5.33 × 105), 5.08 × 105 (3.898 × 105) and 1.78 × 105 (1.43 × 105), respectively, for young, elder and total persons.
Let us estimate the critical number of susceptible persons Sth from equation (19). For R0 = 6.828, we have Sth = 6.532 × 106. Hence, for the São Paulo State, isolating 38.07 million (85.35%) or above persons is necessary to avoid the persistence of epidemics. The number of young persons is 0.27 million less than the threshold number of isolated persons to guarantee eradication of CoViD-19. For R0 = 4.097, we have Sth = 10.87 × 106, implying to isolate 33.73 million (75.63%) or above persons, which is 4.07 million less than young persons.
3.1.2 Fitting the additional mortality rates
We estimate taking into account confirmed deaths from February 26 (t = 0) to April 5 (t = 39), hence n = 40 observations, as we did in the previous estimation.
To estimate the mortality rates αy and αo, we fix the previously estimated transmission rates βy and βo for R0 = 6.828 and 4.097, and evaluate the system of equations (2), (3) and (4), with initial conditions given by equation (7), to calculate
by varying αy and αo. We chose the mortality rates minimizing the sum of differences. Fixing the previously estimated transmission rates βy = 0.75 and βo = 0.8775 (both days−1), we estimate additional mortality rates αy and αo by varying αy and αo. From the fact that lethality among young persons is lower than elder persons, we let αy = 0.2αo [3], and estimate only one variable αo. The estimated rates are αy = 0.0052 and αo = 0.026 (days−1). For lower R0, fixing βy = 0.45 and βo = 0.5265 (days−1), we estimated αy = 0.08 and αo = 0.4(days−1). This is called the first estimation method.
The first estimation method used only one information: the risk of death is higher among elder than young persons (we used αy = 0.2αo). However, the lethality among hospitalized elder persons is around 10% [3]. Combining both findings, we assume that the deaths for young and elder persons are, respectively, 2% and 10% of accumulated cases when Ωy and Ωo approach plateaus (see Figures 4(a) and 5(a) below). This is called the second estimation method, which takes into account a second information besides the one used in the first estimation method. In this procedure, the estimated rates are αy = 0.0018 and αo = 0.009 (days−1) for both R0.
Figure 3 shows the estimated curves of Π = Πy +Πo provided by both methods of estimation and observed data for R0 = 6.828 (a) and 4.097 (b). The second method of estimation fits very badly in the interval of estimation from t = 0 to 39.
The estimated curves of Π provided by both methods of estimation and observed data from t = 0 until 39. The first estimation method provided for R0 = 6.828, αy = 0.0052 and αo = 0.026 (a), and for R0 = 4.097, αy = 0.08 and αo = 0.4 (b). Both figures show also the second method of estimation, which fits very badly in the interval of estimation, providing same αy = 0.0009 and αo = 0.009 (all in days−1).
For higher R0 = 6.828, Figure 4(a) shows the estimated curves of accumulated number of severe CoViD-19 (Ωy, Ωo, and Ω = Ωy +Ωo), from equation (13). At t = 120 days, Ω approached an asymptote (or plateau), which can be understood as the time when the first wave of epidemics ends. The curves Ωy,Ωo, and Ω reach values at t = 200, respectively, 1.511 × 106,0.426 × 106, and 1.937 × 106. Figure 4(b) shows the estimated curves of the accumulated number of CoViD-19 deaths (Πy,Πo, and Π = Πy +Πo), from equation (14). The values of Πy, Πo, and Π are, at t = 200, for the first method of estimation, respectively, 0.747 × 105 (4.9%), 1.137 × 105 (26.7%) and 1.883 × 105 (9.72%), and for the second method of estimation, respectively, 2.67 × 104 (1.77%), 4.767 × 104 (11.19%) and 7.437 × 104 (3.84%). The percentage between parentheses is the ratio Π/Ω.
For R0 = 6.828, the estimated curves of accumulated number of severe CoViD-19 (Ωy, Ωo and Ω = Ωy +Ωo) (a), from top to bottom for total, young and elder persons, and accumulated number of CoViD-19 deaths (Πy, Πo and Π = Πy + Πo) (b), with continuous (αo = 0.026) and dashed (αo = 0.009) curves. At t = 120 days, Ω and Π approached asymptotes (or plateaus).
From Figure 4, the first estimation method for additional mortality rates provided around 2.5-times more than the second method. Especially, the first method estimated 26.7% against 11.2% of deaths in the elder subpopulation, while for the young subpopulation, 5% against 2%. Hence, the second estimation method is more credible.
For lower R0 = 4.097, Figure 5(a) shows the estimated curves of accumulated number of severe CoViD-19 (Ωy, Ωo, and Ω = Ωy +Ωo), from equation (13). At t = 160 days, Ω approached an asymptote (or plateau), which can be understood as the time when the first wave of epidemics ends. The curves Ωy, Ωo, and Ω reach values at t = 200, respectively, 1.484 × 106,0.422 × 106, and 1.906 × 106. Figure 5(b) shows the estimated curves of the accumulated number of CoViD-19 deaths (Πy, Πo, and Π = Πy + Πo), from equation (14). The values of Πy, Πo, and Π are at t = 200, for the first method of estimation, respectively, 6.594×105 (44.4%), 3.583×105 (84.9%) and 1.018 × 106 (53.41%), and for the second method of estimation, respectively, 2.62 × 104 (1.76%), 4.718 × 104 (11.18%) and 7.338 × 104 (3.85%). The percentage between parentheses is the ratio Π/Ω.
From Figure 5, the first estimation method for additional mortality rates provided much higher values than the second method, for instance, 25-times more among young subpopulation. Especially, the first method estimated 85% against 11.2% of deaths in the elder subpopulation. Hence, the second estimation method is more credible. Notice that, from Figures 4 and 5, at the end of the first epidemic wave, we have quite the same number of accumulated CoViD-19 cases for R0 = 6.828 and 4.097 (difference at most 1.8%) despite big differences in accumulated deaths.
For R0 = 4.097, the estimated curves of accumulated number of severe CoViD-19 (Ωy, Ωo and Ω = Ωy +Ωo) (a), from top to bottom for total, young and elder persons, and accumulated number of CoViD-19 deaths (Πy, Πo and Π = Πy +Πo) (b), with continuous (αo = 0.4) and dashed (αo = 0.009) curves. At t = 160 days, Ω and Π approached plateaus.
At t = 0, the numbers of susceptible persons, Sy, So and S = Sy + So are, respectively, 3.77762 × 107,0.68238 × 107 and 4.46 × 107, which diminish due to infection at lower values at t = 200 (figure not shown, see Figure 8(a) below to observe sigmoid shape). For higher R0 = 6.828, the numbers of susceptible persons left behind are 2.346 × 105 (0.62%), 0.294 × 104 (0.043%) and 2.376 × 105 (0.53%), for young, elder and total persons, respectively. For lower R0 = 4.097, the numbers of susceptible persons left behind are 8.823 × 105 (2.34%), 6.957 × 104 (1.02%) and 9.518×105 (2.13%), for young, elder and total persons, respectively. The percentage between parentheses is the ratio S(200)/S(0). For young and total persons, there are 4-times more susceptible persons available to be infected in the second wave when R0 decreases, while for elder persons, 24-times. Hence, the second wave will infect much more elder persons.
For higher R0 = 6.828, the numbers of immune persons (Iy, Io, and I) increase from zero to, respectively, 3.755 × 107 (99.4%), 0.674 × 107 (98.8%) and 4.429 × 107 (99.3%), for young, elder and total persons at t = 200 (figure not shown, see Figure 8(b) below). For lower R0 = 4.097, the numbers of immune persons at t = 200 are 3.689 × 107 (97.7%), 0.668 × 107 (97.9%) and 4.357 × 107 (97.7%), for young, elder and total persons, respectively. The percentage between parentheses is the ratio I/S(0). The value of R0 does not matter in the number of immune persons.
Remembering that human population is varying due to the additional mortality (fatality) of severe CoViD-19, we have, at t = 0, N0y = 3.780 × 107, N0o = 0.680 × 107 and N0 = N0y + N0o = 4.460 × 107, and at t = 200, Ny = 3.774 × 107 (0.16%), No = 0.667 × 107 (1.91%) and N = 4.441 × 107 (0.43%) for the first estimation method, and Ny = 3.779 × 107 (0.026%), No = 0.674 × 107 (0.88%) and N = 4.453 × 107 (0.16%) for the second estimation method. The percentage of deaths is the ratio (N0j − Nj) /N0j.
3.1.3 Estimation of the proportion of isolated persons
From Figures 4 and 5, as R0 decreases, the estimation of the additional mortality rates based on a few initially collected data becomes more inaccurate, providing so much deaths in the population. Hence, we choose βy = 0.75 and βo = 0.8775 (higher R0 = 6.828) and lower estimates for the mortality rates, αy = 0.0018 and αo = 0.009 (all in days−1), to estimate isolation (described by proportions ky and ko) of susceptible persons as control mechanism. Isolation was introduced at t = 27 (March 24), and we will estimate taking into account the confirmed cases from t = 27 to 55 (April 21), hence n = 29 observations.
Here, we estimate the control variables by varying ky and ko. The system of equations (2), (3), and (4), with boundary conditions given by equations (8) and (9), is evaluated, and we calculate the sum of square differences
where t1 = 27 and t29 = 55. We assume k = ky = ko and varied k = 0, 0.4, 0.6, 0.7, and 0.8, from which the better estimated value is k = 0.5 [1]. In this section, the isolation was initiated at t = 27, and release is not initiated.
Figure 6(a) shows the estimated curve of D2 = D2y + D2o and the observed data, plus four curves with lower and higher values of k than that estimated k = 0, 0.4, 0.6, 0.7 and 0.8. Figure 6(b) shows the curves of D2 for 6 different values of k, extended from t = 0 until 250. For k = 0, 0.4, 0.5, 0.6, 0.7 and 0.8, the peaks are, respectively, 6.889 × 105, 3.117 × 105, 2.25 × 105, 1.423 × 105,0.68 × 105 and 0.137 × 105, which decrease as k increases, and displace to right, to t = 68, 82, 88, 99, 118 and 165.
The estimated curve (for k = 0.5) of D2 and observed data, plus five curves with lower and higher values of k than that estimated, k = 0, 0.4, 0.6, 0.7 and 0.8 (a), and curves of D2 extended from t = 0 until 250 (b), from top to bottom as k increases.
Figure 7 shows the curves of Ωy, Ωo, and Ω = Ωy + Ωo (a) and Πy, Πo, and Π = Πy +Πo (b) without (k = 0) and with isolation (k = 0.5). For k = 0.5, the curves Ωy, Ωo, and Ω attain values at t = 250, respectively, 7.296×105,2.09×105 and 9.382×105, and Πy, Πo and Π attain, at t = 250, respectively, 1.29 × 104 (1.77%), 2.334 × 104 (11.17%) and 3.623 × 104 (3.86%). The percentage between parentheses is the ratio Π/Ω.
Figure 8 shows the curves of Sy, So, and S = Sy + So (a) and Iy, Io, and I = Iy + Io (b) without (k = 0) and with isolation (k = 0.5). The numbers of susceptible persons, Sy, So, and S = Sy + So are, respectively, at t = 250 days, 9.948 × 105 (2.63%), 7.942 × 104 (1.17%) and 1.075 × 106 (2.41%), for k = 0.5. The percentage between parentheses is the ratio S(250)/S(0). The number of immune persons (Iy, Io, and I) increases from zero to, respectively, 1.811 × 107 (47.9%), 0.329 × 107 (48.4%) and 2.141 × 107 (48%), for k = 0.5. The percentage between parentheses is the ratio I/S(0), and S(0) can be found in the foregoing section.
From foregoing section, we transport values for k = 0: for Ωy, Ωo, and Ω = Ωy + Ωo, 1.511 × 106,0.426 × 106 and 1.937 × 106; for Πy, Πo and Π, 2.67 × 104 (1.77%), 4.767 × 104 (11.19%) and 7.437×104 (3.84%), ratio is Π/Ω; Sy, So, and S = Sy +So (left behind at t = 200), 2.346 × 105 (0.62%), 0.294 × 104 (0.043%) and 2.376 × 105 (0.53%), ratio is S(200)/S(0); and for Iy, Io, and I,3.755 × 107 (99.4%), 0.674 × 107 (98.8%) and 4.429 × 107 (99.3%), ratio is I/S(0). At t = 0, for Sy, So, and S, we have 3.77762 × 107, 0.68238 × 107 and 4.46 × 107. From Figures 6–8, the peak of severe CoViD-19 cases decreases from 688.9 thousand to 225 thousand and is delayed in 20 days, from 68 (May 4) to 88 (May 24).
The curves of Ωy, Ωo and Ω = Ωy + Ωo (a) and Πy, Πo and Π = Πy +Πo (b) without (k = 0, continuous curves) and with isolation proportion (k = 0.5, dashed curves), from top to bottom for total, young and elder persons.
As k increases to 0.6, 0.7, and 0.8, the peaks decrease to 142.3, 68, and 13.7 (thousand), and are delayed in 31, 50, and 97 (days), occurring at 99 (June 4), 118 (June 23) and 165 (August 9). Severe CoViD-19 cases Ωy, Ωo, and Ω = Ωy + Ωo decreased to 51%, 49% and 48%, while Πy, Πo, and Π, to 48%, 49% and 49%. Susceptible persons Sy, So, and S left behind increased to 424%, 2, 701% and 453%, while Iy, Io, and I, decreased to 48%, 49% and 48%. During the first wave, isolation decreased the number of severe cases and deaths around 50%, as well as immune persons decreased in 50%. However, susceptible persons left behind increased around 5-times. These pictures show that the next second wave will occur earlier and more intense when compared to k = 0.
The curves of Sy, So and S = Sy + So (a) and Iy, Io and I = Iy + Io (b) without (k = 0, continuous curves) and with isolation proportion (k = 0.5, dashed curves), from top to bottom for total, young and elder persons.
Table 3 shows the values of Ω, Π, S, and I at t = 250, for k = 0.6, 0.7, and 0.8. Percentages are calculated with respect to k = 0. As k increases, the numbers of severe CoViD-19 and deaths decrease, which help in the hospital management and care system. However, the remaining susceptible persons increase and less persons are immune, which indicates an anticipation of the second wave. Especially among elder persons, where the number of susceptible persons is so high. Notice that Figure 6 indicates that the number of susceptible persons did not reach the plateau at t = 250 when k = 0.8.
The values of Ω, Π, S and I at t = 250, for k = 0.6, 0.7 and 0.8. y, o and Σ stand for, respectively, young, elder and total persons.
Until now, we estimated the model parameters aiming understanding new coronavirus dy-namics in the São Paulo State. Next, we evaluate the release of isolated persons, which will occur on May 11.
3.2 Epidemiological scenarios considering unique isolation followed by releases
To obtain epidemiological scenarios, we consider the higher R0 and the second method of estimation for additional mortality rates. Additionally, we choose k = 0.5 for the proportion of isolated persons in the São Paulo State. Hence, the fitted values βy = 0.75, βo = 0.8775, αy = 0.0018, αo = 0.009 (days−1), and k = 0.5 are fixed hereafter, unless explicitly cited.
In this section, considering only one isolation occurred on March 24 (t = 27), we study the epidemiological scenarios of release when this isolation will be ended on May 10 (t = 74).1 We stress the fact that the isolation of susceptible persons taken into account by the model discriminates who has or do not have contact with the circulating new coronavirus. In the absence of mass testing, it is expected that among isolated persons there will be who harbor virus. We simulate the beginning of release occurring at t = 72, but also at t = 56 (April 22).2 We consider 3 strategies of release: all released at t = 72 (56) (strategy 1); two releases equally distributed at t = 72 (56) and t = 79 (63) (strategy 2); and three releases equally distributed at t = 72 (56), t = 79 (63) and t = 86 (70) (strategy 3). We also consider two regimes: equal proportions of release for young and elder persons (regime 1), and different proportions (regime 2).
To obtain the scenarios, we solve numerically the system of equations (2), (3), and (4) with initial conditions (t = 0) given by equation (7), the boundary conditions in isolation occurred at t = 27 given by equations (8) and (9), and boundary conditions of releases with first one occurring at t = 72 given by equations (10) and (11).
Aiming comparison with the absence of isolation, when k = 0, the peaks for young and elder persons are 5.024 × 105 and 1.665 × 105, which occur at t = 68. When k = 0.5, the peaks for young and elder persons are 1.667 × 105 (33%) and 5.835 × 104 (35%), which occur at t = 88. Percentage is peak(k = 0.5)/peak(k = 0).
3.2.1 Regime 1 – Equal proportions of releasing young and elder persons
For strategy 1, we have a unique releases with l1j = 1, for strategy 2, two releases with l1j = 0.5 and l2j = 1, and for strategy 3, three releases with l1j = 0.33, l2j = 0.5 and l3j = 1, j = y, o.
Figure 9 shows the curves of D2 without (k = 0) and with (k = 0.5) isolation initiated at t = 27, but released according to strategy 1. For release beginning at t = 56 for young (a) and elder (b) persons, the peaks are, respectively, 4.957 × 105 (99%) and 4.607 × 105 (92%), with peak occurring at t = 82; and for release at t = 72, 1.641 × 105 (99%) and 1.529 × 105 (92%), with peak occurring at t = 92. The release in comparison with only isolation, peaks are anticipated in 6 and delayed at 4 days, respectively, for young and elder persons.
The curves of D2 without isolation (k = 0), with isolation (k = 0.5) initiated at t = 27, but releasing according to strategy 1 beginning at t = 56 (dot and dashed) and 72 (dashed) for young (a) and elder (b) persons. Equal releasing of young and elder persons.
Figure 10 shows the curves of D2 without isolation (k = 0), with isolation (k = 0.5) initiated at t = 27, but released according to strategy 2. For release beginning at t = 56 for young (a) and elder (b) persons, the peaks are, respectively, 4.904 × 105 (98%) and 4.373 × 105 (87%), with peak occurring at t = 84; and for release at t = 72, 1.624 × 105 (98%) and 1.455 × 105 (87%), with peak occurring at t = 94. The release in comparison with only isolation, peaks are anticipated at 4 and delayed in 6 days, respectively, for young and elder persons.
The curves of D2 without isolation (k = 0), with isolation (k = 0.5) initiated at t = 27, but releasing according to strategy 2 beginning at t = 56 (dot and dashed) and 72 (dashed) for young (a) and elder (b) persons. Equal releasing of young and elder persons
Figure 11 shows the curves of D2 without isolation (k = 0), with isolation (k = 0.5) initiated at t = 27, but release according to strategy 3. For release beginning at t = 56 for young (a) and elder (b) persons, the peaks are, respectively, 4.777 × 105 (95%) and 3.97 × 105 (79%), with peak occurring at t = 87; and for release at t = 72, 1.585 × 105 (95%) and 1.33 × 105 (80%), with peak occurring at t = 98. The release in comparison with only isolation, peaks are anticipated at 1 and delayed in 10 days, respectively, for young and elder persons.
The curves of D2 without isolation (k = 0), with isolation (k = 0:5) initiated at t = 27, but releasing according to strategy 3 beginning at t = 56 (dot and dashed) and 72 (dashed) for young (a) and elder (b) persons. Equal releasing of young and elder persons.
Table 4 shows the values of Ω, Π, S, and I at t = 250, for strategies 1, 2, and 3 for release occurring at t = 72. Percentages are calculated with respect to k = 0. When young and elder persons are released simultaneously, the numbers of severe CoViD-19 and deaths maintain practically unchanged, as well as immune persons. However, the remaining susceptible persons increase, which indicates that the second wave harm more elder persons. There is very tiny difference among the three strategies.
Regime 1: the values of Ω, Π, S and I at t = 250, for strategies 1, 2 and 3 for releasing occurring at t = 72. Percentages are calculated with respect to k = 0.
Table 4 presents what happens at the end of the first wave, showing that there is not a significant difference when release is done at t = 56 or 72, either for the three strategies. However, during epidemics, release later (t = 72) is better, as well as the third strategy.
3.2.2 Regime 2 – Different proportions of releasing young and elder persons
Regime 2 deals with releasing young persons, but maintaining elder persons for a while and releasing all them after young persons. This regime presents quite similar figures observed in regime 1. Hence, we illustrate only one case.
Figure 12 shows the curves of D2 without isolation (k = 0), with isolation (k = 0.5) initiated at t = 27, but release young persons according to strategy 3 beginning at t = 56 (dot and dashed) and 72 (dashed) for young (a) and elder (b) persons. All elder persons are released 21 days after the beginning of releasing young persons, that is, at t = 77 and 95. For release beginning at t = 56 for young (a) and elder (b) persons, the peaks are, respectively, 4.67 × 105 (93%) and 3.867 × 105 (77%), which occur at t = 87 and 98; and release beginning at t = 72 for young (a) and elder (b) persons, 1.539 × 105 (92%) and 1.176 × 105 (71%), which occur at t = 92 and 109. The release in comparison with only isolation, peaks are anticipated in 1 and delayed in 10 days (t = 77), respectively, for young and elder persons, and delayed in 4 and 21 days (t = 99) respectively, for young and elder persons. Comparing with Figure 11, the epidemic is slightly mild, but elder persons are infected lately, which could be a desirable effect.
The curves of D2 without isolation (k = 0), with isolation (k = 0.5) initiated at t = 27, but releasing young persons according to strategy 3 beginning at t = 56 (dot and dashed) and 72 (dashed) for young (a) and elder (b) persons. All elder persons are released 21 days after the beginning of release among young persons, that is, at t = 77 and 93.
Now, we assess epidemiological scenarios when elder persons are maintained isolated until end of the first wave t = 250, but young persons are released according to three strategies, that is, for strategy 1, we have a unique l1y = 1, for strategy 2, l1y = 0.5 and l2y = 1, and for strategy 3, l1y = 0.33, l2y = 0.5 and l3y = 1. Elder persons remain isolated, hence lio = 0, for all i.
Figure 13 shows the curves of D2 without isolation (k = 0), with isolation (k = 0.5) initiated at t = 27, but release according to strategy 1. For release beginning at t = 56 for young (a) and elder (b) persons, the peaks are, respectively, 4.81 × 105 and 4.493 × 105, which occur at t = 83 and 92; and for release at t = 72 for young (a) and elder (b) persons, 7.714 × 104 and 6.623 × 104, which occur at t = 83 and 89. There is minor increasing in epidemic among elder persons (peaks are anticipated) due to the increased number of susceptible young persons.
The curves of D2 without isolation (k = 0), with isolation (k = 0.5) initiated at t = 27, but releasing according to strategy 1 beginning at t = 56 (dot and dashed) and 72 (dashed) for young (a) and elder (b) persons. Releasing only young persons.
Figure 14 shows the curves of D2 without isolation (k = 0), with isolation (k = 0.5) initiated at t = 27, but released according to strategy 2. For release beginning at t = 56 for young (a) and elder (b) persons, the peaks are, respectively, 4.764 × 105 and 4.261 × 105, which occur at t = 85 and 95; and released at t = 72 for young (a) and elder (b) persons, 7.532 × 104 and 6.3514.81 × 104, which occur at t = 84 and 90. There is minor increasings in epidemic among elder persons (peaks are anticipated) due to the increased number of susceptible young persons.
The curves of D2 without isolation (k = 0), with isolation (k = 0.5) initiated at t = 27, but releasing according to strategy 2 beginning at t = 56 (dot and dashed) and 72 (dashed) for young (a) and elder (b) persons. Releasing only young persons.
Figure 15 shows the curves of D2 without isolation (k = 0), with isolation (k = 0.5) initiated at t = 27, but release according to strategy 3. For release beginning at t = 56 for young (a) and elder (b) persons, the peaks are, respectively, 4.652 × 105 and 3.865 × 105, which occur at t = 87 and 98; and for release at t = 72 for young (a) and elder (b) persons, 7.276 × 104 and 6.158 × 104, which occur at t = 86 and 89. There is minor increasings in the epidemic among elder persons (peaks are anticipated) due to the increased number of susceptible young persons released.
The curves of D2 without isolation (k = 0), with isolation (k = 0:5) initiated at t = 27, but releasing according to strategy 3 beginning at t = 56 (dot and dashed) and 72 (dashed) for young (a) and elder (b) persons. Releasing only young persons.
Table 5 shows the values of Ω, Π, S, and I at t = 250, for strategies 1, 2, and 3 for release occurring at t = 72. Percentages are calculated with respect to k = 0. When young persons are released while elder persons are maintained isolated, for young classes, the numbers of severe CoViD-19 and deaths maintain practically unchanged, but there are slightly decreasing in immune and increasing in susceptible persons. However, for elder persons there is huge decreasing in the numbers of severe CoViD-19 and deaths, and overall decreasing in severe cases and deaths. There is a tiny difference among the three strategies.
Regime 2: the values of Ω, Π, S and I at t = 250, for strategies 1, 2 and 3 for releasing occurring at t = 72. Percentages are calculated with respect to k = 0.
Comparing Tables 4 and 5, regime 2 shows an advantage with respect to regime 1 by decreasing the number of severe cases and deaths, and, additionally, delaying the peak of epidemics among elder persons. Instead of maintaining isolated indefinitely, we release all elder persons at 4 different times. The first (t = 90) is just after the peak of epidemics without release, the second (t = 120) when occurs inflexion, the third (t = 130) and the fourth (t = 140) when epidemics are ending.
Figure 16 shows the curves of D2 for elder persons without isolation (k = 0), with isolation (k = 0.5) initiated at t = 27, but release of young persons according to strategy 3 beginning at t = 56 (dot and dashed) and 72 (dashed). Disregarding the beginning of release, all elder persons are released at t = 90 (a), 120 (b), 130 (c), and 140 (d).
The curves of D2 for elder persons without isolation (k = 0), with isolation (k = 0.5) initiated at t = 27, but releasing of young persons according to strategy 3 beginning at t = 56 (dot and dashed) and 72 (dashed). Disregarding the beginning releasing time, all elder persons are released at t = 90 (a), 120 (b), 130 (c) and 140 (d).
As elder persons are released lately, they are protected, and the peak decreases and moves rightwardly with respect to young persons, making more easier the efforts of hospitals to face the pandemic outbreaks. However, the second epidemics wave may occur earlier.
Figure 17 illustrates the pulse release of persons by regime 1 for young persons (a) and regime 2 (b), taking into account strategy 3. Curves without (k = 0) and with (k = 0.5) isolation are included to observe difference in the dynamics.
Illustration of pulse release of persons by regime 1 for young persons (a) and regime
2 (b). We show the strategy 3 including curves without (k = 0) and with (k = 0:5) isolation.
4 Discussion
The system of equations (2), (3), and (4) is simulated to provide epidemiological scenarios using parameters estimated from data collected in the São Paulo State, daily released by the Ministry of Health of Brazil. However, these scenarios are more reliable if based on credible values assigned to the model parameters (see discussion in [21]). Here, when we used lower estimated transmission rates, R0 = 4.097, to estimate lethality rates against the data, we obtained αy = 0.08 and αo = 0.4(days−1). Using these estimated fatality rates and lower R0, at the end of the first wave of epidemics, we obtained an incredible 85% of deaths among elder persons.
The substantial difference in the number of deaths relies on the fact that fatality rates are per-capita. In another word, if the number of diseased persons is high, the small fatality rate results in an elevate number of deaths in comparison with the small amount of population. However, the transmission rate is not influenced by the phase of epidemic, due to the fact that this rate multiplies susceptible and infectious individuals.
We assumed that elder persons are under higher risk of being infected by the new coronavirus than young persons, which can be withdrawn if it does not match with actual evidences. Our model did not take into account re-infection (or, maybe relapse), neither if the asymptomatic persons will ultimately manifest symptoms. These unanswered questions do not affect the dynamics during the first wave of epidemics (see [18] for a description of parameters influencing short-term epidemics, and other parameters driving long-term).
Considering that 50% of persons were isolated at the simulation time t = 27 (March 27), but they will be released on May 8 (t = 72), we draw some epidemiological scenarios. We analyzed three release strategies and concluded that there is little difference among them, but the best scenario was obtained if we release young persons maintaining elder persons isolated. In this case, the number of deaths due to CoViD-19 was reduced by half.
In all scenarios of release with k = 0.5, we observed the collapse of the health care system. However, if we increase the proportion of isolated persons above 80%, then the health care system can treat well all patients during the first wave. Let us evaluate the occupancy of beds in hospitals and ICU, using equations (16) and (17). We let h1 = h1y = h1o = 1, h = hy = ho = 0.3, ζ1 = 1/14 (2 weeks of hospital care), and ζ2 = 1/21 (3 weeks of ICU care) [3].
Figure 18 shows the number of beds occupied in hospitals (a) and ICU (b) for k = 0.5 and k = 0.7 (dashed curves), and beds occupied in hospitals (c) and ICU (d) for k = 0.8. For k = 0.5, the peaks of occupied beds in hospitals and ICU (between parentheses) are 1.488 × 105 (8.226 × 104), 4.86 × 104 (2.216 × 104) and 1.895 × 105 (1.043 × 105), occurring at t = 90 (92), 88 (90) and 90 (92), respectively for young, elder and total persons. For k = 0.7, the peaks of occupied beds in hospitals and ICU (between parentheses) are, respectively, for young, elder, and total persons, 4.646×104 (2.728×104), 1.29×104 (0.738×104) and 5.931×104 (3.462×104), occurring at t = 121 (125), 118 (121) and 120 (124). For k = 0.8, the peaks of occupied beds in hospitals and ICU (between parentheses) are 9, 495 (5, 879), 2, 674 (1, 597) and 12, 160 (7, 471), occurring at t = 170 (176), 165 (169) and 169 (174), respectively for young, elder and total persons.
The number of beds occupied in hospitals (a) and ICU (b) for k = 0.5 and k = 0.7 (dashed curves), and beds occupied in hospitals (c) and ICU (d) for k = 0:8.
For k = 0.5, the peaks of occupied beds by all persons in hospital and ICU are 189, 500 and 104, 300, occurring at t = 90 (May 26) and 92 (May 28). For k = 0.7, the peaks of occupied beds by all persons in hospital and ICU are 59, 310 and 34, 620, occurring at t = 120 (June 25) and 124 (June 29). For k = 0.8, the peaks of occupied beds by all persons in hospital and ICU are 12, 160 and 7, 471, occurring at t = 169 (August 13) and 174 (August 18). In this case, From Figure 18, the peak of severe CoViD-19 cases is 13, 700 occurring at t = 165 (August 9).
As we have pointed above, when isolation is done lately, among isolated persons, there will be persons harboring virus. Just after the beginning of isolation (t = τis+), the number of persons in each class, for k = 0.5, is
with Q1y = Q1o = 0 and I = 8.8376 × 103. To estimate the circulation of virus among isolated persons, we simulate the system of equations (2), (3), and (4) taking as initial conditions the number of persons in each class, at t = τis+: Sy = Qy = 1.8858 × 107, So = Qo = 3.4056 × 106, Qy = Qo = 0, and for all other variables, half of the corresponding values at t = τis+. Figure 19 shows the curves of D2y, D2o, and D2 = D2y + D2o for 
, with R0 = 1.37 (a), and 
, with R0 = 0.68 (b), from t = 27 to 74, during the period of isolation. At t = 72, the numbers of severe CoViD-19 cases are, 1, 597, 635 and 2, 231, respectively, for young, elder, and total persons for R0 = 1.37; and 294, 129 and 423, for R0 = 0.68. For k = 0.5 and R0 = 1.37, this 0.8% additional CoViD-19 cases compared to the peak is negligible, however for k = 0.8, the 18% additional CoViD-19 cases are not negligible. Nevertheless, when these isolated infectious individuals are released, they contribute to increase the velocity of propagation.
The curves of D2y, D2o and D2 = D2y + D2o for 
(a) and 
(b) from t = 27 to 74, during the period of isolation.
The virus should be circulating among isolated persons through restricted contact occurring in the household and/or neighborhood. If symptomatic persons are transferred to hospital, then removed from isolation, it is expected that asymptomatic persons are predominantly spreading the new coronavirus. If we can assume that they are releasing a low amount of virus, it is expected that more asymptomatic cases will arise than severe CoViD-19. Increasing immune persons, they decrease susceptible persons, hence the effective reproduction number decreases, which is a herd immunity phenomenon.
Let us discriminate the circulation of the new coronavirus in a community. Figure 20 shows all persons harboring this virus (Ej, Aj, D1j, Q2j, and D2j), j = y, o, for young (a) and elder (b) persons corresponding to Figure 6(b) with k = 0.5. There is little difference: exposed (E) and pre-diseased (D1) persons are relatively higher among young persons.
The curves of all persons harboring this virus (Ej, Aj, D1j, Q2j and D2j), j = y; o, for young (a) and elder (b) persons.
In Figure 21 we show the ratio hidden:apparent based on Figure 20. Who harbor the new coronavirus as exposed or who do not manifest symptoms are classified in the hidden category, and in the apparent category, we include who manifests symptoms. Hence, the ratio is calculated as (E + A + D1):(Q2 + D2). At t = 0, the ratio is 10 : 1 for young and elder persons due to initial conditions.
The curves of ratio hidden:apparent for young, elder and total persons.
Comparing Figures 20 and 21, as the epidemics evolves, the ratio increases quickly in the beginning, reaches a plateau during the increasing phase, and decreases during the declining phase, finally reaching another plateau after the ending phase of the first wave. In the first plateau, the ratios are 18 : 1, 26 : 1, and 24 : 1 for, respectively, elder, young, and total persons. The second plateau (5 : 1) is reached when the first wave of epidemics is ending. Therefore, during epidemics, there are much more hidden than apparent persons, which makes any control mechanisms hard if mass testing could not be implemented. However, the estimation of the ratio between hidden and apparent cases can be helpful in designing mass testing to isolated asymptomatic persons.
In this paper, epidemiological scenarios were obtained aiming the evaluation of isolation and subsequent release. However, from Figure 6(a), we observed a lowering in severe CoViD-19 data with respect to the estimated curve. Hence, let us assume that at t = 45, 18 days after the beginning of isolation, protection actions were adopted by persons (use of alcohol gel and face mask) [6]. Figure 22(a) shows three curves of D2 = D2y + D2o and observed data, where 
, 0.5βy = 0.375, 0.45βy = 0.3375 and 0.4βy = 0.3. The better estimated is 
. Figure 22(b) shows the curves of D2 for 4 different decreasing values of transmission rates extended from t = 0 until 250. We stress the fact that this evidence must be confronted with more data.
Three curves of D2 = D2y + D2o and observed data, where 
, 0.45βy = 0.3375 and 0.4βy = 0.3 (a). The better estimated is 
. Curves of D2 for 4 different decreasing values of transmission rates extended from t = 0 until 250 (b)
The question of using a face mask and constant hygiene (washing hands with alcohol and gel and protection of the mouth, nose and eyes) to avoid infection is left to future work.
5 Conclusion
We formulated a mathematical model considering two subpopulations comprised by young and elder persons to study CoViD-19 in São Paulo State, Brazil. The model considered pulses in isolation and release. Briefly, the first CoViD-19 case was observed on February 26 (t = 0), isolation were initiated on March 24 (t = 27), and the release will be implemented on May 8 (t = 72).
We estimated transmission rates, additional mortality (lethality) rates, and proportions iso-lated from data collected in the São Paulo State. To estimate transmission rates and additional mortality (lethality) rates, we used data collected from February 26 to April 5 (t = 39), and data collected from March 24 to April 21 (t = 55), to estimate the proportion isolated. We estimated high values for the basic reproduction number, R0 = 6.828, which is encountered among airborne viruses (see [14] for an estimation of R0 = 6.71 for rubella infection). Currently, some authors found SARS-CoV-2 aerosol in some areas of the hospital environment, indicating the possibility of infection through air [22].
Using the estimated parameters, the system of equations (2), (3), and (4) were simulated to providing epidemiological scenarios when release will be implemented after initial isolation. Using the estimated proportion k = 0.5 of isolated persons, we observed that in all three strategies of releasing for regime 1 (equal release of young and elder persons), the severe CoViD-19 cases approach those without isolation, but the peaks are delayed. However, for regime 2 (releasing young but maintaining elder persons), deaths due to CoViD-19 are reduced by half among elder persons, and by 30% in all persons. Nevertheless, simulation showed that an isolation of 80% is desirable aiming reduction in the severe CoViD-19 cases.
In our model, we did not allow severe CoViD-19 cases to transmit the infection, as well as persons with mild symptoms also. This isolation assumption may have an indirect effect on the transmission of disease by two hypotheses: (1) if a virulent strain is causing severe cases, then their isolation must decrease its fitness, and (2) if severe cases release more amount of virus, then their isolation may decrease the abundance of virus in the environment. As a consequent, asymptomatic and mild cases of CoViD-19 can prevail as epidemics evolves. However, immediate consequent will be infection with more virulent strain and more absorption of the virus by heath care workers, which is a preeminent reason to provide them with extremely secure equipment.
Finally, severe CoViD-19 data collected in the São Paulo State indicates that there is another lowering in those cases besides the diminishing resulted from isolation. We hypothesize that this decreasing could be due to protection actions, but more data are needed to confirm or deny this hypothesis.
Data Availability
The data that support the findings of this study are openly available in Ministry of Health (Brazil) at https://covid.saude.gov.br.
A Trivial equilibrium and its stability
By the fact that N is varying, the system is non-autonomous non-linear differential equations. To obtain autonomous system of equations we let kj = lij = 0, j = y, o, and use fractions of individuals in each compartment, defined by, with j = y and o,
resulting in
using equation (5) for N. Hence, equations (2), (3) and (4) in terms of fractions become, for susceptible persons,
for infected persons,
and for immune persons
where λ is the force of infection given by equation (1) re-written as
and
which is autonomous system of equations. We remember that all classes vary with time, however their fractions attain steady state (the sum of derivatives of all classes is zero). This system of equations is not easy to determine non-trivial (endemic) equilibrium point P*. Hence, we restrict our analysis with respect to trivial (disease free) equilibrium point.
The trivial or disease free equilibrium P0 is given by
for j = y and o, where
with 
.
Due to 17 equations, we do not deal with characteristic equation corresponding to Jacobian matrix evaluated at P0, but we apply the next generation matrix theory [4].
The next generation matrix, evaluated at the trivial equilibrium P0, is obtained considering the vector of variables x = (ey,ay,d1y,eo,ao,d1o). We apply method proposed in [15] and proved in [16]. There are control mechanisms (isolation), hence we obtain the reduced reproduction number Rr by isolation.
In order to obtain the reduced reproduction number, diagonal matrix V is considered. Hence, the vectors f and v are
and
where the superscript T stands for the transposition of a matrix, from which we obtain the matrices F and V (see [4]) evaluated at the trivial equilibrium P0, which were omitted. The next generation matrix FV−1 is
and the characteristic equation corresponding to FV−1 is
where the basic reproduction number R0 is
and R0y and R0o are the basic partial reproduction numbers defined by
Instead of using the spectral radius 
, we apply procedure in [15] (the sum of coefficients of characteristic equation), resulting in a threshold R0. Hence, the trivial equilibrium point P0 is locally asymptotically stable (LAS) if R0 < 1.
In order to obtain the fractions of susceptible individuals, M must be the simplest (matrix with least number of non-zeros). Hence, the vectors f and v are
where superscript T stands for the transposition of a matrix, from which we obtain the matrices F and V evaluated at the trivial equilibrium P0, which were omitted. The next generation matrix FV−1 is
and the characteristic equation corresponding to FV−1 is
The spectral radius is 
given by equation (A.8). Hence, the trivial equilibrium point P0 is LAS if ρ < 1.
Both procedures resulted in the same threshold, hence, according to [20], the inverse of the reduced reproduction number R0 given by equation (A.8) is a function of the fraction of susceptible individuals at endemic equilibrium s* through
where 
(see [18] [20]). For this reason, the effective reproduction number Ref [17], which varies with time, can not be defined neither by Ref = R0 (sy + so), nor Ref = R0ysy + R0oso. The function f (ϰ) is determined by calculating the coordinates of the non-trivial equilibrium point P*. For instance, for dengue transmission model, 
, where 
and 
are the fractions at equilibrium of, respectively, humans and mosquitoes [18]. For tuberculosis model considering drug-sensitive and resistant strains, there is not f (ϰ), but s* is solution of a second degree polynomial [20].
From equation (A.10), let us assume that 
. Then, we can define the approximated effective reproduction number Ref as
which depends on time, and when attains steady state (Ref = 1), we have s* = 1/R0.
The basic partial reproduction number 
(or 
) is the secondary cases produced by one case of asymptomatic individual (or pre-diseased individual) in a completely susceptible young persons without control; and the partial basic reproduction number 
(or 
) is the secondary cases produced by one case of asymptomatic individual (or pre-diseased individual) in a completely susceptible elder persons without control. If all parameters are equal, and ψ = 1, then
where 
and 
are the basic partial reproduction numbers due to asymptomatic and pre-diseased persons.
The global stability follows method proposed in [9]. Let the vector of variables be x = (ey,ay,d1y,eo,ao,d1o), vectors f and v, by equations (A.5) and (A.6), and matrices F and V evaluated from f and v at trivial equilibrium P0 (omitted here). Vector g, constructed as
results in
and gT ≥ 0 if 
and αy = αo = 0.
Let vl = (z1,z2,z3,z4,z5,z6) be the left eigenvector satisfying vlV−1 F = ρvl, where 
,
and
This vector is
and Lyapunov function L, constructed as L = vlV−1xT, is
always, and
only if 
and αy = αo = 0.
Hence, the method proposed in [9] is valid only for αy = αo = 0, in which case P0 is globally stable if 
and 
.
Footnotes
↵a hyunyang{at}ime.unicamp.br,
↵b luispedro_jr{at}hotmail.com,
↵c ffmcastro{at}gmail.com,
↵d arianacy{at}gmail.com
↵1 Initially, programmed to end on April 6, further extended to April 22, and postponed to May 10.
↵2 Simulations were done on April 21. We anticipated the release to May 8 (May 10 is Mothers’ day).
