In this work we give a delay diferential equation , the retarded logistic equation, as a mathematical model for the global transmission of COVID-19. Th is model accounts for asymptomatic carriers, pre-symptomatic or latent transmission as well as contact tracing and quarantine of suspected cases. We find that the equation admits varied classes of solutions including self-burnout, progression to herd immunity and multiple states in between. We use the term “partial herd immunity” to refer to these states, where the disease ends at an infection fraction which is not negligible but is significantly lower than the conventional herd immunity threshold. We believe that the spread of COVID-19 in every localized area can be explained by one of our solution classes.
Retarded logistic equation
• Applied computing – mathematics and statistics
Three kinds of models to study COVID-19 are currently in vogue – lumped parameter or compartmental models (ordinary differential equation), agent-based models and stochastic differential equation models. The first option affords maximum conceptual clarity at the expense of some simplifying assumptions †Presenting author, Corresponding author. ORCID : 0000-0003-2502-2268 (homogeneous mixing etc). The second option affords maximum potential versatility at the cost of huge computational complexity and variability in the network structure. The third option combines features of the previous two – whether the features being synergized are the positive or the negative ones depends to a large extent on the modeler.
In this work we use delay differential equations (DDE) to
propose a simple, single-variable, lumped parameter model for the
spread of Coronavirus. Jahedi and Yorke [
We measure time t in days and use as our basic variable y(t) which is the cumulative number of corona cases, including active cases, recovered cases and deaths, in the region of interest. The following “word-equation” summarizes the approach : Rate of emergence Interaction rate of = of new cases each existing case Probability of Number of transmission existing cases persons/day. For a compartmental model, one must average over the professor, the grocer and all the other un-quarantined cases to generate an effective per-case interaction rate q0.
Every interaction of course does not result in a transmission – there is a probability strictly less than unity that the virus jumps from the infected person to the person whom s/he is interacting with. This probability has two components. The first component is that the healthy person must be susceptible to begin with. While we ignore intrinsic insusceptibles, there will be people who have recovered from the disease and are therefore not susceptible again. In this Article, we assume that one bout of infection brings permanent immunity. The assumption is valid so long as the immunity period exceeds the total epidemic duration. Till date, there is little credible evidence for re-infection [5]–[7]; contrarily, a very recent and thorough study [8] based on monitoring of huge patient cohort has found significant evidence of long-lasting and effective antibodies. If N be the initial number of susceptible people (recall that y is the case count), then the probability that a random person is a recovered case is approximately y/N and the probability that s/he is susceptible is (approximately) 1−y/N. This expression is approximate because the true number of recovered cases at any time is less than y; the error however is small since the recovery period is much shorter than the overall course of the epidemic. Note that 1−y/N is a logistic term, and a herd immunity effect.
Given susceptibility, the next probability is that the virus actually does jump from the un-quarantined case to the susceptible person. This probability depends on the level of precaution such as face covering or mask, handwashing and disinfection being adopted by the case as well as the susceptible person. For a compartmental model, the probability must be averaged over all the un-quarantined cases. If this average probability is P0, then q0(1−y/N)P0 gives the per-case spreading rate. Since q0 and P0 are both dependent on public health measures, and are both difficult to measure independently, we can club those two together into a single parameter which we call m0.
So far we have accounted for the rate at which each cases spreads the disease; now we have to count the number of cases out of quarantine. Let us start with an asymptomatic carrier, who remains in open society throughout. S/he typically transmits the disease for 7 days, which is called the infection period. Then, new healthy people can be only be infected by those asymptomatic cases who have fallen sick within the last 7 days, and not those who have fallen sick earlier. The number of such people is the number of asymptomatic sick people today minus the number of those 7 days earlier. Mathematically, let μ1 (between 0 and 1) denote the fraction of asymptomatic carriers and τ1 the asymptomatic infection period. Then, the number of asymptomatic transmitters today is μ1(y(t)−y(t−τ1)). Here we can see the emergence of the delay term.
The remaining fraction 1−μ1 of cases are symptomatic. Let τ2 be the latency period during which these cases remain transmissible prior to displaying symptoms. It is assumed that they isolate themselves thereafter. Assumption is also made that the incubation period is equal to the latency period. Finally, the
Due to complexity of the equation (2), analytical solution using perturbation theory etc has not been attempted in this case. Instead we have used numerical integration to obtain the solutions of (2). Before giving the solutions however, we present the calculation of the reproduction number R. To find R at any state of evolution of the disease, we first treat y in the logistic term to be constant, and then carry out the steps described in Ref. [9]. This yields the expression
R = m (1 − 0 y ) ( 1 + μ3 − 2 μ1 μ3 τ + μ μ τ )
2 1 3 1
N 2
The ease of calculating R with respect to the ordinary differential equation based models [10] is noteworthy.
Solution classes of logistic DDE (2) are now demonstrated. The numerical integration routine used is second order Runge Kutta with a time step of 1/1000 day. As the testbed for the simulations, we consider a Notional City having N=300000, μ1=0.8, (maximum . (3) contact tracing drive conducted by public health department is taken into account. Assumption is made that this drive is instantaneous and proceeds in forward direction starting from freshly arriving symptomatic cases. The contact trace captures patients who were exposed to the new case τ2 days ago, as well as patients who were exposed immediately before the new case manifested symptoms. The average duration for which these secondary patients have remained at large is τ2/2, be they symptomatic or asymptomatic. The assumption of instantaneous contact tracing, which decreases the average time that contacttraced cases spend out of quarantine, opposes the error arising from the assumption of zero non-transmissible incubation period, which increases the average time for which the contact-traced cases transmit before quarantine. These two effects are assumed here to cancel. Let μ3 (between 0 and 1) denote the fraction of all cases who escape from contact tracing drives – the complementary fraction 1−μ3 get caught. Thus, we have three classes of un-quarantined cases : (a) 1−μ3 are contact-traced cases who remain in society for a time τ2/2, (b) μ3 (1−μ1) are untraced symptomatic cases who go into isolation only after time τ2, and (c) μ3μ1 are undetected asymptomatic cases who transmit for the entire infection period τ1. Arguments similar to those of the previous paragraph yield the total number of un-quarantined cases as n = (1 − μ3 ) ( y − y (t − τ 2 / 2) ) + ( (1 − μ1 ) μ3 ) ( y − y (t − τ2 ) ) + μ1 μ3 ( y − y (t − τ1 ) ) The preceding arguments now yield the mathematical form of (0) as dy dt = m
0 1 − y y (t ) − (1 − μ3 ) y (t − τ2 / 2) − N (1 − μ1 ) μ3 y (t − τ2 ) − μ1 μ3 y (t − τ1 ) (2) which is the retarded logistic equation.
. (1) value as per our knowledge [11]–[13]), τ1=7 days and τ2=3 days [14]. The initial condition needs to be a function having the length of the maximum delay involved in the problem, which is seven days; we take this function to be zero cases to start with and constant increase of 100 cases/ day for a week.
Notional City A has m0=0.23 and μ3=1/2, which describes a hard lockdown [15] accompanied by good contact tracing. R0 (i.e. (3) evaluated at y=0) is 0.886. The epidemic ends with a negligible fraction of infected people, as shown below. This and the next five plots are three-way – each plot shows y as blue line, its derivative y as green line and the weekly increments in cases, or epidemiological curve, as a grey bar chart. These last have been reduced by a factor of 7 to ensure clarity of presentation. We report the rates on the left hand side y-axis and the cumulative cases on the right hand side y-axis.
This is exactly what has happened in New Zealand – that il fortunatissimo per verita has indeed quashed the epidemic completely with the final case count being a negligible fraction of its total (tiny and sparsely distributed) population.
The parameter values for Notional City B are the same as those for A except that μ3=0.75; a greater fraction of cases escape the contact tracing drive. R0 is 1.16, and R becomes 1 at y=40500 cases.
The outbreak enters exponential regime right after being released. As y increases, R gradually reduces so the growth slows down until it peaks when the case count is about 39,000 [compare with the value of 40,500 when R=1 as per (3)]. Thereafter, the disease progresses to extinction in time. The overall progression is very long but one hopes that the relatively small size of the peak can prevent overstressing of medical care facilities and thus avoid unnecessary deaths. Delhi and Mumbai in India and Los Angeles in USA are in all probability cities of this type since the disease there spiraled out of control despite hard lockdowns being imposed at an early stage.
City B also enables us to explain partial herd immunity. Even though the initial conditions were unfavourable for containment of the epidemic, herd immunity started activating as the disease proliferated. A stable zone (R<1) was entered when only 13.5 percent of the total susceptible population was infected, and a similar percentage again got infected before the epidemic ended. Thus, herd immunity worked in synergy with nonpharmaceutical interventions to stop the epidemic at only 26 percent infection level, which is significantly less than the conventional 70-90 percent threshold [16]. This is what we call partial herd immunity. Our findings are in agreement with and act as an explanation for what has been obtained by Britton et. al. [17] and Peterson et. al. [18].
We now consider Notional City C which differs from City B in that m0=0.5; lockdown is replaced by a much more permissive state. R0 is above 2.5; 1,80,000 infections are required to bring it below unity.
Figure 3 : City C goes to herd immunity – total not partial. The symbol ‘k’ denotes thousand and ‘L’ hundred thousand.
Need one mention that this is a public health disaster. Notional City D combines features of B and C. This city begins with m0=0.5 like City C but reduces to m0=0.23 like City B when the case count reaches 40,000 (the R=1 threshold for B’s parameters).
Figure 4 : As the input, so the output – D’s response combines features of B and C. The symbol ‘k’ denotes thousand and ‘L’ hundred thousand.
We can see a case count as well as a total duration intermediate to B and C; the epidemic is over in 70 days but the peak rate of 12,920 cases/day is still very high and likely to load hospital facilities beyond their carrying capacity.
The Cities E and F demonstrate the issues faced in reopening. In both these cities, the parameters and case trajectory are identical to those of City A for the first 80 days. Then, E and F reopen on the 80th day by increasing m0 from 0.23 to 0.5, and simultaneously decreasing μ3 i.e. deploying a more effective contact tracing program which had been built up during the lockdown. The post-reopening μ3’s for E and F are 0.1 and 0.2 respectively.
The difference between Cities E and F is dramatic. Mathematically, R remained less than unity throughout in E; its value after reopening was 0.985. We can see that the case rate decreases monotonically all the time. In F, the post-reopening R became 1.22 and sent the trajectory haywire. In practice however, the incipient increase in case rate after the 80th day acts as an advance warning of what has happened – the reopening steps should be reversed if it is at all possible to do so while satisfying economic and other external constraints.
In this Article we have presented a new mathematical model for COVID-19 which is simple and elegant in structure but can generate a variety of realistic solution classes. We hope that our work may be of use to mathematicians and data scientists who are trying to understand the spread of the disease in a quantitative manner. The public health implications of these results are being reserved for another study.