Abstract
COVID-19 is caused by a hitherto nonexistent pathogen, hence the immune response to the disease is currently unknown. Studies conducted over the past few weeks have found that the antibody titre levels in the blood plasma of infected patients decrease over time, as is common for acute viral infections. Fully documented reinfection cases from Hong Kong, India, Belgium and USA, as well as credible to anecdotal evidence of second-time cases from other countries, bring into sharp focus the question of what profile the epidemic trajectories may take if immunity were really to be temporary in a significant fraction of the population. Here we use mathematical modeling to answer this question, constructing a novel delay differential equation model which is tailored to accommodate different kinds of immune response. We consider two immune responses here : (a) where a recovered case becomes completely susceptible after a given time interval following infection and (b) where a first-time recovered case becomes susceptible to a lower virulence infection after a given time interval following recovery, and becomes permanently immunized by a second infection. We find possible solutions exhibiting large number of waves of disease in the first situation and two to three waves in the second situation. Interestingly however, these multiple wave solutions are manifest only for some intermediate values of the reproduction number R, which is governed by public health intervention measures. For sufficiently low as well as sufficiently high R, we find conventional single-wave solutions despite the short-lived immunity. Our results cast insight into the potential spreading dynamics of the disease and might also be useful for analysing the spread after a vaccine is invented, and mass vaccination programs initiated.
Introduction
Over the last month, one fully documented case of COVID-19 reinfection has been discovered in Hong Kong [1], two such cases in Greater Noida, India [2,3], one more in Belgium [4] and a fifth in Seattle, USA [5]. In general, the second infections appear to be milder than the first. These discoveries intensify the relevance of questions such as the duration and degree of immunity conferred by one bout of infection. Due to the rapidly evolving nature of the disease, the unknowns in this area outnumber the knowns by multiple orders of magnitude. A past study [6] on benign human coronaviruses (not SARS-CoV, MERS-CoV or SARS-CoV-2!) has found that immunity against infection lasts for a few months and reinfection is common from one year onwards. In a recent study [7], a research group from Mount Sinai Hospital in New York, USA found that among a cohort of almost 20,000 patients, all but one demonstrated significant levels of antibody titre in blood plasma at three months following the original infection. These results were corroborated by a smaller study conducted at Max Hospital, New Delhi, India [8].
An observational cohort study [9] on 34 patients conducted at the University of Washington at Seattle found that over a 3 to 5 month period, antibody concentrations decreased over time, in a manner consistent with the immune responses to acute infections by other viruses including influenza, SARS and MERS. In these infections, the initial decrease in titres is followed by a plateau; whether this is true of SARS-CoV-2 is obviously unknown as yet. This study also reports that, while severe infection caused the initial peak of antibody concentration to be significantly higher than mild or asymptomatic infection, these higher-concentrated antibodies also decayed much faster so that at the end of three months, there was little noticeable difference between the titre levels of severe and asymptomatic patients. A study conducted at Doha, Qatar [10] has reported that out of more than 1,30,000 patients infected with COVID-19, 243 persons reported a positive swab at least 45 days after the original positive test; 54 out of these 243 had “strong or good evidence for reinfection”. All the reinfections were asymptomatic or had minimal or mild symptoms on the second bout; however the initial symptom profile of these patients had been mild or asymptomatic as well. A recent study emerging from Thailand [11] reports upto 20 percent post-infection seroconversion failure rate in mild or asymptomatic patients – it does not document any reinfection cases. Some potential cases of reinfection have also been reported in media for some time before the first documented case – in these instances the evidence is not fully credible [12-14].
Two recent systematic reviews deal with the antibody response [15] and the cellular immune response [16] to the virus. The former contains amplifications of the kind of results we discussed above; the latter basically states that the response of T-cells to SARS-CoV-2 is currently unknown. Finally, a recent study [17] proposes that the immune response to the virus will determine the manner of transition of COVID-19 to endemicity; in the absence of known facts, this work is predominantly speculative at this time.
Although the reinfection cases so far are isolated, it is natural to wonder what the epidemiological consequences might be if the immunity duration indeed turns out to be finite for a significant fraction of the population. The only approach which can allow inroads into this question is mathematical modeling. Giordano et. al. [18] and Bjornstad et. al. [19] account for the possibility of reinfection, with the latter finding an oscillatory approach towards the eventual endemic equilibrium. Kosinski [20] however finds multiple waves of COVID-19 if the immunity threshold is finite. Before commencing our own analysis of this situation, let us clarify that we are currently treating large-scale temporary immunity as a hypothetical, possibly worst case scenario, regarding the validity of which we do not yet have sufficient available data.
Methods
We start from a delay differential equation (DDE) model constructed earlier by our group [21,22] which accounts for many realistic features associated with COVID-19 transmission. The dependent variable y(t) denotes the cumulative number of corona cases in the region of interest (typically a town, neighbourhood, village or other area with good interaction among the inhabitants) as a function of time. The parameters in the model are :
m0 : the per-case spreading rate which accounts for factors such as degree of mobility (i.e. lockdown/unlock), extent of mask use, extent of handwashing and other public health measures
τ1 : the asymptomatic transmission period which we take to be 7 days throughout
τ2 : the pre-symptomatic transmission period which we take to be 3 days throughout
μ1 : a number between 0 and 1, it denotes the fraction of patients who are asymptomatic
μ3 : a number between 0 and 1, it denotes the fraction of patients (both symptomatic and asymptomatic) who are NOT detected in contact tracing drives and quarantined
N : the total number of susceptible people in the region at the start of the epidemic
In Refs. [21,22], we had assumed that one bout of infection renders everlasting immunity (practically, immunity lasting longer than the epidemic’s progression). Now, we account for two different kinds of immune response, as follows.
Simple response : In this scenario, one bout of infection renders a person insusceptible to fresh infection for a time τ0 days (τ0 much greater than τ1 and τ2) following recovery, after which s/he again turns susceptible to the disease in its original form
Complex response : In this scenario, adapted from Ref. [17], a person is initially susceptible to the current, highly virulent form of the disease. The first bout of infection renders him/her completely insusceptible for τ0 days following recovery, after which s/he turns susceptible to a lower virulence form of the disease. If infected the second time, then s/he becomes permanently insusceptible to further infection
As is customary in lumped parameter or compartmental models, the value of τ0 used in the model must be an average over the entire population. The experiences of the first few confirmed reinfections may well be indicative of a general trend favouring the second immune response to the first; in the absence of data, we shall work with both assumptions.
It can be shown that with the simple immune response, the dynamics of the disease is governed by the DDE
The derivation of this equation is given in the Supplementary Material. With the complex immune response, we need two dependent variables, y(t) the cumulative number of cases of the high virulence form and z(t), the cumulative number of cases of the lower virulence form. We let μ1a and μ1b (presumably greater than μ1a) be the asymptomatic fractions for the two forms and keep all other parameters identical for both forms. Then, the dynamics of y and z are governed by the following coupled system of DDEs
whose derivation is again given in the Supplement.
For both models, we present the solutions in a Notional City having an initial susceptible population of N=3,00,000 and 80 percent asymptomatic carriers i.e. μ1=0.8 (or μ1a=0.8 if appropriate) in the current (high virulence) form of the disease. Six such cities corresponding to six fundamental solution classes have been demonstrated in Refs. [21,22] under the assumption of permanent immunity. We let μ1b=0.95 for the lower virulence form and take the immunity threshold τ0 to be 200 days. We solve the equations numerically in Matlab using 2nd order Runge Kutta method with step size 0.001 day. The seeding function we take is zero cases for the first 193 days followed by linear growth of cases at 100 cases/day for the next seven days (for complex immune response, this growth refers to the high virulence cases – the lower virulence cases remain zero). We take t=0 to be the 194th day of the seeding period.
One issue needs special mention – in a numerical simulation, the case rate will never be identically zero but will be something like 0.001 cases/day (or machine epsilon cases/day). This can pose a serious problem in a situation where there are potential second waves of epidemic. To circumvent this issue, we have arranged for manual termination of the run if the case rate becomes low enough. Defining the number of active cases at time t to be y(t)− y(t− 14), we stop the run if there is less than one active case for 14 consecutive days. While the number 14 (twice) may be somewhat arbitrary, the criterion of a low enough active case count for a long enough period is a very reasonable indicator of the true end of the outbreak. We run all simulations either upto t=1400 days or until they terminate, whichever is earlier.
Results
All results in this Section are based on an assumed immunity threshold of 200 days.
Simple immune response
We first consider the simple immune response, modelled by (1). In Refs. [21,22] where immunity was taken as permanent, Notional City A had m0=0.23 and μ3=0.5, which resulted in the epidemic’s being driven monotonically to containment in 120 days. With the simple immune response, the time-trace of the disease is shown in Fig. 1 below. Here and henceforth, we show three things in the same plot : y(t) as a blue line, its derivative as a green line and the “epi-curve” or weekly increments in cases scaled down by a factor of 7, in grey bars.
City A drives the epidemic to extinction in time.
The response remains the same as it was earlier – since the epidemic was contained starting with everyone susceptible, it does not matter whether a fraction have since become immune and how long they remain so.
Next, we present City C which has μ3=3/4 and m0=0.5 all the time. (We use 0.23 for a “low” value of m0 since we obtained it from data fits in Refs. [21,22]; the chosen “high” value of m0=0.5 generates 90 percent infection level at the end of the outbreak.) With permanent immunity, City C went up to this high infection level in just 50 days. With simple immune response, the trajectory is shown in Fig. 2 below.
City C reaches herd immunity well before the 200-day mark. ‘k’ denotes thousand and ‘L’ hundred thousand.
There is no change in behaviour, on account of the short duration of the epidemic in City C. City B has m0=0.23 like City A but μ3= 3/4 like City C. With permanent immunity, City B crawls up to about 26 percent infection level over a period of 220 days. Here however, the 220-day run becomes longer than the assumed immunity threshold of 200 days, and the result, now labelled as City B1, is shown in Fig. 3.
City B1 (we shall use B2 to denote the same city with complex immune response) has multiple waves of COVID-19. ‘k’ denotes thousand and ‘L’ hundred thousand. We have stopped the simulation at 1400 days – the waves persist after that time as well.
We can see wave after wave of outbreaks here. While we do not believe that the disease can really run unmitigated for four years, we have chosen a long simulation runtime to show the perfect periodicity of the case trajectories. The case count at the end of the run is greater than the city’s initial susceptible population, so at least some people have been infected at least twice.
Cities D and E of Refs. [21,22] were similar to C and A respectively and again they do not show any change when the new immune response is incorporated. City F, which demonstrates reopening, is more interesting. We start it off with the parameter values of B, on a path to multiple waves. 150 days into the outbreak, it reopens completely raising m0 to 0.5, aiming to infect its entire population before immunity runs out. The result is shown in Fig. 4.
City F is bold and “beats the virus to the finish line”. ‘k’ denotes thousand and ‘L’ hundred thousand.
Howsoever controversial this strategy might appear, it works. The reopening generates an immediate second wave but it does succeed in making the epidemic vanish completely at 210 days and 2,46,000 infections. This strategy however has an unhappy variant – City G, instead of reopening at one shot, increases m0 linearly from 0.23 to 0.50 over the interval from 150 to 250 days. The resulting infection profile is shown in Fig. 5.
City G1 is F gone wrong, and the error carries a price. ‘k’ denotes thousand and ‘L’ hundred thousand.
There are three separate waves of infections, each more severe than the first, and the cumulative case count exceeds 150 percent of the total population.
Complex immune response
We now consider the complex immune response, modelled by (2). Numerical work confirms our expectation that Cities A, C and F will show no difference from the simple immune response case since the epidemics there terminate before the immunity duration lapses. The first City we focus on is B, now renamed to B2. For this plot the colour legend will be 
in blue, y in green, z in red, ż in magenta and the epi-curves in grey for y and cyan for z. The results are in Fig. 6.
City B with complex immune response. ‘k’ denotes thousand and ‘L’ hundred thousand. W/7 denotes the weekly increments in cases scaled down by 7 as previously, while HVF and LVF refer to the high and lower virulence forms of the disease respectively.
Once more we can see the multiple waves but this time they are progressively attenuated so that the total fraction of the high virulence infection remains less than 60 percent.
An equally dramatic difference between the two immune responses occurs with City G, Fig. 7 below.
City G with complex immune response. In the right third of the plot, the cyan bars obscure the smaller grey bars; the latter are enveloped by the green curve which is still visible. ‘k’ denotes thousand and ‘L’ hundred thousand. W/7 denotes the weekly increments in cases scaled down by 7 as previously, while HVF and LVF refer to the high and lower virulence forms of the disease respectively.
The massive third wave of G1 now occurs almost entirely in the lower virulence form.
Discussion and conclusion
The difference between the cities A, B and C lies in their reproduction number R, which is proportional to m0 if other parameters are held constant [21,22]. City A has the starting value R0=0.886, City B has R0=1.16 and City C has R0=2.5. With permanent immunity, a low but greater-than-unity R0 leads to a long epidemic duration with a lower caseload while high R0 leads to a shorter duration with higher caseload. The latter very likely results in overwhelming of healthcare facilities and causes unnecessary deaths. With temporary immunity however, we are seeing the possibility of multiple waves of disease if the threshold τ0 becomes less than the evolution time with permanent immunity (see the Supplementary Material for some pedagogical examples which motivate these results). In this example, we took τ0 to be 200 days, which is towards the longer end of the free evolution duration. Hence we see the wave solutions in a small region of parameter space – we find that with the simple immune response, for all parameters except m0 being held to their City B values, a containment (City A) solution for low m0 gives way to a multiple wave (City B1) solution as m0 is increased above 0.20, while that cedes to a one-shot logistic-like (City C) solution at m0=0.25 and higher. A shorter immunity duration would have caused multiple waves over a larger range of m0. With the complex immune response, the total numbers of cases are bounded by 3,00,000 infections of both y and z. Hence, the infinitely periodic waves are no longer possible – either they attenuate in size or the epidemic terminates. For example, if m0 in City B2 is raised slightly to 0.24, then the disease stops at 550 days with two waves and total 1,70,000 high virulence and 35,000 low virulence infections (this m0 in City B1 causes infinitely many waves). In both Cities B and G, the time interval between successive waves appears to be in the range 1.5τ0 to 2τ0.
We note that the second waves seen after relaxing public health interventions are fundamentally different from the second and subsequent waves being exhibited under the aegis of temporary immunity. The former is caused by more susceptible people getting exposed and infected, as in the second wave of City F, while the latter is caused by recovered cases turning re-susceptible and thus increasing the size of the susceptible pool, like the waves in City B1. The waves in City G1 are of both types. The second wave is primarily the result of the reopening starting on the 150th day while the third wave is the result of the first and second wave people having turned re-susceptible. This is why in City G2, where previously infected people suffer a “different” disease, the second wave occurs primarily in the high virulence form while the third occurs in the lower virulence form. The phenomenon we are currently seeing in European countries and Israel is overwhelmingly the result of weakening of public health intervention measures and not of immunity having run out. We are not sure as to what factors governed the successive waves of the deadly 1918-9 influenza pandemic, and how that disease eventually ended.
We now present a comparison between the results obtained here and in the prior Literature. Some References [18,19,23,24,25] appear to have missed the multiple wave solutions possible with temporary immunity. Reference [19] finds an oscillatory approach to an endemic equilibrium, while Ref. [23] finds steady oscillations about such an equilibrium – both features being absent in conventional S-E-I-R models. Ref. [20] gets the multiple waves like we do. But it finds these waves whatever the value of R0, while we obtain diverse possible responses depending on the value of R0. We believe that the present results are more realistic, as explained in the Supplement. Finally, we have never seen an S-E-I-R framework being used to model a situation such as the complex immune response where a person can be infected exactly twice but not more.
In conclusion, our work is speculative as of now since it is predicated on an assumption whose validity is not known at the current time. Given the extreme variety of the solutions we find, our analysis probably adds to the body of unknowns and uncertainty about this new disease. Still, it is better to be prepared in advance for every scenario instead of being blindsided by it in real time. A further application of our study will arise when a vaccine is invented and mass vaccination initiated – then also, the immunized population can be modelled by our process and the trajectories of the disease can be estimated. We shall elaborate this aspect further after a vaccine is actually released for mass use.
Data Availability
There is NO data referred to in the manuscript.
Funding statement
We have NO funding to report for this study.
Conflict of interest statement
We have NO conflict of interest.
Author contribution statement
Both of us contributed equally to the manuscript.
SUPPLEMENTARY MATERIAL
Derivation of the model equations (1) and (2)
The philosophy behind the delay differential equation model has been proposed in Refs. [S1,S2] of this Supplement (Refs. [21,22] of the Article proper). We start by summarizing the derivation presented in these References. This is necessary because we shall adopt this derivation only to obtain (1) and (2) of the Article proper.
Recapitulation of the prior works [S1,S2]
The transmission of COVID-19 arises as a result of interactions, either direct or via objects, between cases at large in society and healthy people. It is described by the following “word equation” :
which we intend to represent mathematically.
The first term on the above right hand side (RHS) is itself the product of two quantities : the interaction rate q0 of each at large case with other people, and the probability P0 that an interaction between a case and a susceptible person actually results in a transmission. q0 is determined by the level of social mobility (lockdown/unlock) and by the degree of physical separation being practised; P0 depends on whether either or both parties are wearing masks and whether they are washing hands, not touching their face etc. Both q0 and P0 depend directly on public health interventions, so we can combine them into a single term m0 which represents this aspect of the disease.
The second term on the RHS of (S0) is the probability that the person with whom the at large interacts i.e. a random person in the region, is actually susceptible. For now we consider the case where immunity is permanent. In this case, the probability of a random person’s being insusceptible is the total number of recovered cases divided by the total number of people in the region. The former is approximately y, and the latter is N. The recovery count is actually a little less than y since cases take a finite time to recover. But if the recovery duration is much shorter than the overall epidemic duration, which it is in practice, then we can make the assumption of instantaneous recovery and treat it as y. With this assumption, the probability of a random person’s being insusceptible is y/N and the probability of his/her being susceptible is 1− y/N.
The structure of the third term arises from the following argument : if all cases take a time τ days to recover and remain at large the whole time, then the number of at large active cases now is exactly the number of people who fell sick between now and τ days back – mathematically, this is expressed as y(t)− y(t− τ). For a maximally realistic picture, we partion the cases into three classes : contact traced cases, untraced symptomatic cases and untraced asymptomatic cases. By definition of the model parameters, the first class amount to fraction 1− μ3 of the total cases; if the contact tracing starts from freshly reporting symptomatic cases, then the average duration these cases remain at large is τ2/2 where τ2 is the transmissible latency period. Untraced symptomatic cases account for fraction μ3(1− μ1) of the total cases and they remain at large for the latency period τ2 before manifesting symptoms and going into quarantine. Finally, untraced asymptomatic cases account for fraction μ3μ1 of the total cases and these remain at large for the asymptomatic infection period τ1. Using the y(t)− y(t− τ) argument on each class leads to the mathematical form of the third term on the RHS of (S0) as
which is the retarded logistic equation proposed in Refs. [S1,S2].
For further details regarding the derivation of (S2), including justification of the approximations involved, we must refer to Refs. [S1,S2]. We shall now use the above arguments to derive the model equations (1) and (2) of the Article proper, which form the backbone of the present contribution.
Construction of the model
We briefly recapitulate the two kinds of immune response from the Article proper. Simple immune response : After contracting the infection, a person remains insusceptible for τ0 days following recovery and then becomes susceptible again. Complex immune response : After contracting the infection, a person remains insusceptible for τ0 days, then becomes susceptible to the lower virulence form and after a second infection becomes permanently immune.
We first consider the simple immune response. In the structure (S0), immune response is modelled by the second term on the RHS, so that is the only term which will be different from (S2). The immunity threshold of τ0 means that at any given instant, the people who are immune are all those who have contracted the infection during the last τ0 days, and no one else. Hence, the number of insusceptibles at time t is exactly the number of new infections which have occurred between time t− τ0 and t, which is y(t)− y(t− τ0). In proposing this structure, we have retained the approximation of instantaneous recovery and also have ignored deaths. Since the mortality rate of COVID-19 is fortunately quite low, this second assumption is reasonable as well. So the probability that at time t, a random person is insusceptible is [y− y(t− τ0)]/N, and the probability that s/he is susceptible is its complement,
Using this in (S0) immediately gives
which is (1) of the Article proper.
With complex immune response, we need to write the structure (S0) for both y and z. The public health term m0 remains the same for both. In the y-equation, the second term must denote the probability that a random person is susceptible to the high virulence form. As before, let us count the insusceptibles. Any person who has contracted the high virulence form once is insusceptible to it for all future time, so the number of insusceptibles is y and the probability of susceptibility is 1− y/N. For the third term, we need the total numbers of active and at large y as well as z, which is (S1) written out for both y and z. Thus, we get
which is (2a) of the Article proper.
For the z equation, the per-case spreading rate as well as the number of spreaders remain unchanged; the only thing which changes is the probability that a random person is susceptible to the low virulence form. At time t, the only people eligible for contracting the low virulence form are those who have already contracted the high virulence form once, and that a time τ0 or more ago. This is the cumulative number of high virulence infections at time t− τ0, which is y(t− τ0). Among these eligible people however, the current count z(t) have already had the second infection and so must be excluded from the susceptible pool. Thus, the size of the susceptible pool is y(t− τ0)− z(t) and the susceptibility probability is
leading to the equation
which is (2b) of the Article proper. This completes the model derivation.
Qualitative motivation and explanation of the solutions
The different solution classes of (1) and (2) we saw in the Article proper are quite counter-intuitive so we motivate them here through some pedagogical examples. We had used one such example in an earlier work [S3] to motivate the City A containment solution, which we now recapitulate briefly. This example was the fantasy kingdom of Covidland where the king declares a 28-day “100 percent lockdown” – literally nobody is allowed to step out of his/her house. Because Covidland is fictitious, this does not pose a logistical problem of survival. At the end of the lockdown, every case has either recovered or died at home, and the virus does not exist any longer when normal life is restarted. City A is the real-world version of Covidland since the philosophy is the same – with skilful public health measures, the virus is fed new patients at such a small rate that it itself gets annihilated in time.
The pedagogical examples relevant for temporary immunity are the following. The fictitious city of Covidtown has an initial population of 500 susceptible people, and has no cases to start with. The immunity threshold is 200 days and every case recovers within a day (consistent with the modeling assumption above). The spread of the virus here is governed by the following rule :
On any day if there are 5 or more than 5 susceptible individuals present, then exactly 5 of them (randomly selected where necessary) contract the virus
On any day if there are less than 5 susceptible individuals present then all of them contract the virus
The case trajectory of Covidtown is simple enough to predict. There are total 5 cases on day 1, 10 on day 2, 15 on day 3 and so on, all the way upto 500 on day 100. Then, since the first batch of cases is still well within its immunity period, there are no more susceptible people left and the propagation stops. On day 100, the last batch of cases also clears the virus which then doesn’t exist in Covidtown any longer. Everyone has been infected but the virus is dead and the epidemic over.
The fantasy city of Covidpolis has everything same as Covidtown except that the number 5 is replaced by 2. The case trajectory starts off with 2 cases on the first day, 4 on the second, 6 on the third, 200 on the 100th (unlike the 500 of Covidtown) and so on upto 400 cases on the 200th day. Thus, on day 201, there are 100 people who are yet to contract the infection. But, on this day, the cases of day 1 lose their immunity and also get added to the susceptible pool. So the two new infections of day 201 can be anyone from this pool of 102. Similarly on day 202, the two cases from day 201 get removed from the susceptible pool but the two recoveries from day 2 get added and the pool size does not change. This scenario remains invariant for all time – every day there is a susceptible pool of 102 and every day there are 2 new infections from amongst this pool. Even though the initial growth rate is much slower than in Covidtown, the city of Covidpolis becomes a land of immortal virus.
The above is exactly what we see in the solutions of (1) and (2). A high R like in Cities C and F leads to a fast spreading of the epidemic with high initial case counts but puts the progress of the disease to a stop in time. A lower R (but of course greater than unity) leads to slower spreading with lower initial case counts but the disease keeps on perpetuating itself, as in City B1. These examples help us to qualitatively understand the solutions and also convince us that our results are a more accurate representation of reality than literature items which do not find this complex dependence on R.
References
