Reformulating the susceptible-infectious-removed model in terms of the number of detected cases: well-posedness of the observational model.

Compartmental models are popular in the mathematics of epidemiology for their simplicity and wide range of applications. Although they are typically solved as initial value problems for a system of ordinary differential equations, the observed data are typically akin to a boundary value-type problem: we observe some of the dependent variables at given times, but we do not know the initial conditions. In this paper, we reformulate the classical susceptible-infectious-recovered system in terms of the number of detected positive infected cases at different times to yield what we term the observational model. We then prove the existence and uniqueness of a solution to the boundary value problem associated with the observational model and present a numerical algorithm to approximate the solution. This article is part of the theme issue 'Technical challenges of modelling real-life epidemics and examples of overcoming these'.

Introduction

At the beginning of any pandemic, mathematical modelling of infectious diseases is thrust into the limelight for governments across the globe to gain imperative forecasting pictures of the impending spread of the infectious disease, as demonstrated by the 2009–2010 swine flu pandemic [1-4], the 2015–2016 Zika virus epidemic [5-8] and the most recent COVID-19 pandemic [9-12]. Non-pharmaceutical interventions (NPIs) are interpreted by adjusting model parameters to allow for quick decision-making on the appropriate NPI measure(s) to be implemented. Understanding and predicting levels of infections in the community is the crux of mathematical modelling in epidemiology since it is these infections that drive the epidemic and also the mathematical models. However, getting accurate, reliable data to describe the current number of infections within the community is arguably the largest challenge for data-driven modelling. Not only is lack of data problematic when deriving and using mathematical models, the context and interpretation of the data and the parameter fitting process often comes with difficulties. In particular, there is often not a one-to-one correspondence between data and parameters or compartments in a model, especially when considering daily snapshots of an epidemic, for example.

In this work, we study the susceptible–infectious–removed (SIR) equations often accredited to the work of Kermack & McKendrick in 1927 [13]. In particular, we look to demonstrate how to use the number of detected cases directly in the model formulation and prove the existence and uniqueness of a solution to the resulting formulation. The SIR equations take the following form of a simple temporal dynamic system of ordinary differential equations, supported by non-negative initial conditions where the dot above the notation denotes the time derivative, see figure 1 for a schematic. As is standard for the SIR equations, we let denote the total population being considered, denote the average transmission rate and denote the average removal rate. The initial conditions are chosen so that , which gives rise to for any . Moreover, it is standard to make the assumption that , i.e. there is at least one infected individual at the onset of the infectious disease. We direct readers to the following literature to gain insight into the analysis and uses of the SIR equations and its family of compartmental models [14-17].

Figure 1

Schematic of the SIR compartmental model. We let denote the proportion of the total population who are susceptible to the infectious disease being studied. Susceptible individuals become infectious with the disease to form the subpopulation at a rate , which represents the current average infection rate. The rate is the product of , the average transmission rate, and the probability of meeting an infectious person . Individuals in the subpopulation then are removed at a rate to form the removed subpopulation . The description of the removed compartment depends on the nature of the infectious disease and the interpretation of the data available to be used. Typically, the term removed is interchanged with the term recovered. As is standard for epidemiological models of this nature, denotes the average time between transmissions and denotes the average length of time before removal.

The ‘observational model’, as coined in [18], is a representation of a compartmental model, such as the SIR equations, described in terms of the parameters and compartments that are captured, or to be inferred, by the mathematical interpretation of the data. This observational model gives an intuitive understanding of what parameters can be identified and how the data affects different compartments in the model. Typically, data for infectious diseases come as a daily, or in general, periodic, snapshot of some of the states that individuals experience during the infectious disease. In this study, we are looking to use the detected cases of an infectious disease, which we denote by . In practical terms, we are assuming that our observation is a fraction of the total number of individuals leaving the infectious state over the time of the observation. This situation arises when the individual ceases to be infectious once they have been diagnosed, for example because the individual can avoid contact with susceptible individuals. We allow for the situation that we only observe a fraction of cases since there are situations where not all infectious individuals can be removed, such as false-negatives in testing regimes, bad individual compliance in behaviour or financial/resource restrictions on the amount of testing that can be done. In reality, is likely to be time-dependent since it can characterize aspects of public policy, however, we keep it constant in this study. In order to describe these observations in terms of the model parameters and compartments, we formally describe the removed compartment as those who are no longer mixing in the general population. Hence, for each observation of the number of detected cases , , we see that where represents the time interval of the measurement of the data. This definition of the observed data for the number of detected cases of an infectious disease has been used in the past; see for instance [19].

Using (1.4) as the interpretation of the available data leads to the observational model a second-order nonlinear differential equation equipped with non-local boundary conditions (1.4), see §3 for the detailed derivation. A significant part of this work is devoted to studying the boundary value problem (1.4) and (1.5). In contrast to initial value problems that typically arise when studying mathematical epidemiology, boundary value problems are much more difficult to deal with. For instance, it is well known that a small change in the boundary conditions may lead to significant changes in the behaviour of the solutions [20]. In fact, there is no general theory to the existence and uniqueness of a solution for boundary value problems. In most instances, we have to study the boundary value problem on a case-by-case basis, especially when the problem at hand is nonlinear, and/or the boundary conditions are atypical, which is the case with our problem as the boundary condition (1.4) is given in a non-local integral form. Non-local boundary conditions have been studied in the past in various forms, recently appearing in the study and application of fractional differential equations [21]. In our case, the non-local boundary conditions correspond to fixing the values of certain integrals of the solution over given time intervals, which is similar to how the average value is imposed on certain boundary value problems with Neumann boundary conditions to ensure uniqueness of the solution. Although this form is convenient for the analysis, we note that we can convert the integral conditions to nonlinear, multi-point Robin-type conditions. This, however, is only possible due to the set-up of the observational model and thus the problem at hand, since the differential equation directly involves the data, and thus the boundary condition in some capacity. As far as we are aware, this is the first work to consider a boundary value problem with these specific integral-type boundary conditions. The key message of our study is on the transformation of an epidemiological initial value problem to a boundary value problem due to the appropriate form of the observed data, hence we refer readers to [20-26], and the references therein, for further information on methods and techniques that may be used to analyse certain classes of boundary value problems with different types of boundary conditions.

The main results of this publication can be summarized in the following two theorems:

Theorem 1.1 (Existence).

Let . Given corresponding to ( , , there exists a solution to (.

Theorem 1.2 (Uniqueness).

Let . Given corresponding to ( , , the solution to ( . Furthermore, the solution is smooth.

The structure of the paper is as follows: in §2, we motivate how the approach outlined in this work came about and why it is necessary, in §3, we derive the observational model and pose an equivalent boundary value problem to (1.5), in §4, we prove Theorem 1.1, in §5, we prove theorem 1.2 and in §6, we outline an algorithm for the numerical solution to the observational model with two examples.

Motivation

Since the onset of COVID-19, the disease caused by SARS-CoV-2, in December 2019 and the declaration of a pandemic by the World Health Organization in March 2020, governments and organizations across the globe have implemented unprecedented regional, national and international interventions to contain the spread of COVID-19 and the subsequent damage caused by it. Initially, NPIs were the only tool in the fight to contain the spread of the disease, forcing governments to introduce strict control measures on public mixing, the use of personal protective equipment, improvement of personal hygiene and the use of contact tracing to reduce the possibility of transmission throughout their countries, protecting those most vulnerable as well as maintaining control of the healthcare demand and capacity. The social-economical effects of these necessary decisions are still unknown, and it is speculated that these effects will still be felt for many more years to come [27-30]. The global picture is starting to look less bleak with the introduction of effective vaccines and vaccination policies allowing the relaxation of such strict NPIs, however with the spawn of new variants of COVID-19, such as the Alpha variant (lineage B.1.1.7) and more recently the Delta variant (lineage B.1.617.2), it is clear we are not out of the pandemic yet [31-34]. Indeed, now that ‘normality’ has returned, country borders have opened and tourism is on the increase, we in the UK are amidst an increase in daily positive cases, hospitalizations and deaths, while other European countries are tightening their COVID-19 restrictions or even going back into lockdowns.

In response to the pandemic in the Sussex region of the UK, local authority organizations in Sussex created a consortium that posed modelling questions of strategic operational significance to the Local Health Resilience Partnership covering East Sussex, West Sussex and Brighton & Hove. As a result, leaders from Public Health Intelligence in the local authority organizations approached the University of Sussex to undertake COVID-19 epidemiological modelling that is specific to Sussex as the national level modelling from the Scientific Advisory Group for Emergencies was not applicable at the regional, integrated care level. In response to their request, we developed a data-driven SIR-type mathematical model, which in practice has been used to answer public health questions based around healthcare demand and capacity, and mortuary capacity, such as the impact of second waves, future lockdowns and vaccination supply [18]. In our early parameter estimation attempts, it became clear to us that the data provided were not compatible with the standard SIR-type formulation, in particular, the data did not provide an accurate estimation of all the initial conditions of the compartments of the model. Healthcare data typically come in one of three types: new daily data (such as admissions), current daily data (such as beds occupied) or cumulative daily data (such as total number of deaths). Due to the methods and constraints on data collection, data typically come as new daily data, i.e. a change in the state of an individual. For example, being admitted into hospital is a change of state from being infectious in the community to being in a hospital bed, whereby the infectious state and hospital bed state are compartments of a mathematical model, while the admitted component is the change between compartments. In this set-up, the admissions data are proportional to the infectious compartment, but not directly as it is an integral, so identifying the initial conditions of the compartments from the data becomes a real challenge, see [18] for more details.

Although the remit of this paper does not include parameter estimation results, we nevertheless describe these in order to demonstrate the issue of the estimation of initial conditions described above. Given that we have a system of nonlinear differential equations, the parameter estimation problem can be stated as a nonlinear least-squares (NLS) optimization problem [35]. Let describe the observed data, describe the set of parameters we are estimating, describe the set of solutions to (1.1) and (1.2) with describing the differential equations, and denote the function that maps the solutions to the mathematical interpretation of the data. Let describe the likelihood function associated with the probability distribution assumption of the noise between and . Then, the NLS optimization problem can be stated as where denotes the admissible parameter space. The NLS problem easily demonstrates the issue of the initial conditions faced: that even although parameters can be estimated regardless of the initial conditions, they may not be representative of real life. This impacts the use of epidemiological models, that they may not have any accurate realistic forecasting capabilities. Moreover, it can easily be reasoned that is highly non-convex with respect to the parameters and in practice will therefore likely converge to solutions close to the initial guess, which may not be the global solution [36,37].

The standard approach in the literature to navigate the issue of estimating initial conditions is to restrict the parameters, for example, by fixing the basic reproductive number, and estimate when the epidemic began using the data provided, see for instance [38]. Fortunately for modellers in the UK, with regards to COVID-19, the disease was brought to us from outside our borders and so identifying ‘patient zero’ and tracing their connections becomes easier. However, with epidemiological data in particular, what is often not considered is that there is stochastic perturbation in the realization of the actual epidemic. The spread of an infectious disease is a stochastic process, and so the collection of data from an epidemic is just one realization of that process. Noise in the data presents itself mathematically as uncertainty in the parameters, while stochastic perturbation in the epidemic presents itself mathematically as uncertainty in both the initial conditions and the parameters. Variability due to stochastic perturbation is high at the onset of an epidemic, especially in a spatially complex and heterogeneous population not typically accounted for in modelling, which makes estimating the first infection and calibrating to data difficult. In [18], we took an alternative approach and estimated the initial conditions from the data at the onset of the first national lockdown, since this was the first available bit of data. In this instance, the initial conditions became part of the parameters we looked to estimate, however this proved challenging due to the nonlinearity of the NLS problem and the differences in parameter sensitivity between initial conditions and system parameters (such as and ). This led us onto the problem we are trying to solve currently: Can we incorporate the data into the model formulation in a more natural way that overcomes the issue of accurate estimation of initial conditions? We believe that by rewriting the SIR equations as the observational model with boundary conditions, we can gain accurate estimates of parameters while bypassing issues with initial conditions. However, before approaching the parameter identifiability and estimation problems, we first need to demonstrate the non-trivial well-posedness of the observational model.

Derivation of the observational model

Before we derive the observational model, it is important we provide some standard analytical results of (1.1) and (1.2). It is a standard result in dynamical systems and phase space analysis that for positive initial conditions, we have and for . In particular, steady-state analysis tells us that where . Due to (1.1), is monotonically decreasing, which implies that In addition, since the right-hand sides of (1.1) and (1.2) are smooth functions, their solutions are infinitely differentiable (see, for instance, [39]), which allows us to make the necessary steps to derive the observational model. Turning our attention to the data (1.4), for the analysis, we assume that, for any , That is to say that we always observe a strictly positive, finite number of cases since we are not considering any observation error or model error in this study, i.e. our observations are exact. Indeed, in this situation, the case implies that there no more infectious cases in the population, and in consequence, there will be no more cases in future observations.

As described in [18], since the data (1.4) are described by the infectious compartment , we need to rewrite equations (1.1) and (1.2) to be purely in terms of in order to derive the observational model. To this end, we note three properties of equations (1.1) and (1.2), namely which, by differentiating, immediately gives us Moreover, by differentiating (1.1), we see that and hence, inserting (3.2) and (3.3) into (3.4), the observational model takes the form as seen in (1.5). One notes that (3.4) and thus (1.5) are well defined for due to (3.1). For the purposes of this paper, we continue to transform (1.5) to facilitate the analysis. By denoting and setting , (1.5) is equivalent to while the data (1.4) are equivalent to In the subsequent sections, we will prove the following two theorems, which are respectively equivalent to Theorems 1.1 and 1.2 due to the continuous bijective transformation between and .

Theorem 3.1 (Existence).

Given , corresponding to ( , there exists a solution to (.

Theorem 3.2 (Uniqueness).

Given corresponding to ( , the solution to ( . Furthermore, the solution is smooth.

Remark.

One can observe that as tends to 0, becomes negatively unbounded. This is not impactful on the analysis due to the boundary condition (3.6), equivalently (1.4), always being positive. As we noted above, in the situation that (3.6) is zero, and upon assuming no noise in the data, this implies that the epidemic is over and thus this model is not relevant anymore.

One can go further to simplify (3.5) by taking , where , to get where the prime denotes the derivative with respect to , which leads to (3.6) satisfying

To continue to demonstrate the issues with initial conditions, we can express and in terms of and the parameters. By definition of , we obtain and From the formulation of the observational model (3.5), we can see that a solution is defined only in terms of and . Assuming that is known, due to the definition of , the pair cannot be uniquely identified from the data; in other words, for a given , any pair satisfying will produce exactly the same solution. This implies that in order to simulate future predictions with the original SIR formulation, one would need estimates for , as well as the initial conditions. Indeed, from (3.7), we conclude that it is not possible to identify unless both and are known; in particular, since it is not possible to identify both and from the data, it is not possible to recover the initial conditions without an independent estimate of the under-reporting coefficient . By contrast, an advantage of the formulation of the observational model is that it will allow us to compute the solution , and predict future observations , as long as we can fit the model to identify the value of . Note also that, without an independent estimate on , we can only infer the ratio between the initial number of infectious individuals to the initial number of susceptible individuals since

Existence

The proof of theorem 3.1 is based on the Leray–Schauder fixed point theorem [40]. We shall show that the solutions to (3.5) and (3.6) satisfy a priori bounds, and we will use them to show compactness of a suitable operator.

Lemma 4.1.

Let and . Suppose there exists a solution to ( and . Then, there exists some constant such that , for . Moreover, is only in terms of the parameters and the data.

Proof.

Let and let . Integrating (3.5) over we have Let and let . Then, by integrating with respect to , we have Since the exponential function is positive everywhere, we see that and consequently Hence, we see that there exist bounds for , where and These bounds depend on and through their values at . Since we are interested in uniform bounds for any solution , we need to show that and are bounded in terms of the data and . By exponentiation and integration (4.1) over , we have and hence, by taking (4.2) into account, To show that and are bounded, we proceed by contradiction. Let and , and let A standard application of L’Hôpital’s rule yields where we have fixed in the first limit and fixed in the second limit, and similarly Since (4.4) implies boundedness of , we conclude and cannot be arbitrarily large and in consequence, and are bounded from above. Using a similar argument with the right-hand side bound of (4.4) and taking the limits and , we can conclude that and are also bounded from below.

Therefore, the bounds (4.2) and (4.3) are independent of the solution , and depend only on the parameters and the data. ▪

We can now proceed to show that is also uniformly bounded in .

Lemma 4.2.

Let and . Suppose there exists a solution to ( and . Then, there exists some constant such that , for . Moreover, is only in terms of the parameters and the data.

Using the notation introduced in the proof of lemma 4.1, we integrate (3.5) over , where , to obtain We can now use the uniform bounds on , from lemma 4.1, together with the bounds on and that we obtained in the proof of lemma 4.1, to bound the right-hand side of (4.5). ▪

We can proceed using equations (3.5) and (3.6) and the bounds on and to obtain uniform bounds on the second derivative, which in turns allows us to obtain uniform bounds for the third derivative by using them in (3.5). Indeed, we can bootstrap this argument to obtain uniform bounds up to the th derivative, for finite . We summarize this in the following corollary.

Corollary 4.3.

Let , with , and let and . Suppose there exists a solution to (3.5) and (3.6) with data and . Then, there exists some constant , depending only on , the data and the parameters, such that for all .

We can now use the uniform bounds on the solutions to study a suitable operator that will allow us to define the solutions of (3.5) and (3.6) as a fixed point. Let where depend on and are uniquely determined under the constraints We note importantly that a fixed point of is a solution to (3.5)–(3.6). Therefore, our goal is to apply the Leray–Schauder fixed point theorem to to show the existence of a solution in . We will first show that is a compact operator.

Lemma 4.4.

Let be defined as in ( is a compact operator in .

Let be a sequence of functions bounded in . We see that since we have is also uniformly bounded, since is uniformly bounded in . Therefore, the Arzela–Ascoli theorem implies that is a compact operator in . We note that to apply the Arzela–Ascoli theorem in suffices to show uniform boundedness of the third derivative, since then we also have equicontinuity of the second derivative; this is analogous to using uniform boundedness of the first derivative to show equicontinuity of the sequence of functions. ▪

With this set-up, in particular that is compact in , we can now finish the proof of theorem 3.1.

Proof of theorem 3.1.

In order to apply the Leray–Schauder fixed point theory, we now need to show that the set is bounded for all . It is easy to see, by linearity of the derivative and the integral, that a fixed point of corresponds to a fixed point of with a change of parameters, namely, by replacing by . Therefore, any function in (4.8) is a solution to (3.5) and (3.6) with a suitable parameter, and as shown in lemma 4.1, is bounded in . Thus the Leray–Schauder fixed point theorem implies the existence of a fixed point for , which, as already mentioned, corresponds to a solution of (3.5)–(3.6). ▪

Uniqueness

Now that we have proven that there exists a solution to (3.5)–(3.6), we now want to demonstrate the uniqueness of a solution to (3.5)–(3.6). We employ a standard contradiction argument for boundary value problems by showing that two different solutions, which satisfy the same boundary conditions, must be the same.

Proof of theorem 3.2.

Let and be two solutions to (3.5) and (3.6) satisfying the same data conditions, that is, for we have Consequently, this means that Let and . By taking the difference of (3.5) for each solution, we have To simplify the presentation, we set and denote Using this notation, (5.2) is equivalent to Let . Integrating (5.3) with respect to gives Moreover, by (5.1) and (5.4), one can see that We now proceed by contradiction to show that . Without loss of generality, assume that . This gives us that by definition and by (5.3). Considering this, (5.5) and the fact that , , we have that there exists such that and for all . In particular, this means that, by integrating with respect to , we have and, setting in (5.4), we have Then, by integrating (5.4) with respect to and using (5.6), we obtain where the last equality is a consequence of Fubini’s theorem and (5.7). Therefore, we have a contradiction. One notices that, using the same approach, we can also conclude that due to (5.5). By integrating (5.4) with respect to , we have that which, by noting that implies that and using (5.1), gives us Standard ODE theory implies that, since , and , , we have for all . ▪

Numerical examples

Now that we have demonstrated existence and uniqueness, we can build a numerical algorithm to approximate the solution to (3.5)–(3.6). We set up the numerical problem as follows: given , , and , we look to find the function that satisfies (3.5)–(3.6). We use the shooting method to numerically approximate the solution to (3.5), see for example [41] for details. We use a superscript 0 to denote the initial approximation to and a superscript to denote the converged solution from the shooting method.

For our first example, we look at the case where the number of detected cases is reducing. We set , and use the data and with , and . For the shooting method, we use an initial guess of , , depicted by the dashed blue line in figure 2a, which results in The resulting initial conditions from the shooting method are and , with the resulting function depicted by the continuous orange line in figure 2a. Even although we can typically gain knowledge of the total population , through for example a population census, as previously mentioned we cannot extract from due to the under-reporting parameter . However, for the purpose of demonstration, we fix and to allow us to demonstrate and , the number of infectious people transformed from and , respectively. Figure 2b depicts the transformed initial guess as the dashed blue line and depicts the transformed converged solution as the continuous orange line.

Figure 2

Demonstration of the numerical approach to approximate the solutions to (3.5) and (3.6) in the case of the number of detected cases is reducing. (a) Comparison of the initial guess for and the converged numerical approximation for . (b) Comparison of the transformed initial guess for and the transformed converged numerical approximation for . (Online version in colour.)

For our second example, we look at the case where the number of detected cases is increasing. We set , and use the data and with , and . For the shooting method, we use an initial guess of and , depicted by the dashed blue line in figure 3a, which results in The resulting initial conditions from the shooting method are and , with the resulting function depicted as the continuous orange line in figure 3a. Again, we fix and to allow us to demonstrate and in figure 3b, whereby the transformed initial guess is depicted by the dashed blue line and the transformed converged solution is depicted by the continuous orange line.

Figure 3

Demonstration of the numerical approach to approximate the solutions to (3.5) and (3.6) in the case of the number of detected cases is increasing. (a) Comparison of the initial guess for and the converged numerical approximation for . (b) Comparison of the transformed initial guess for and the transformed converged numerical approximation for . (Online version in colour.)

Conclusion

In this publication, we have demonstrated how to derive the observational model given data that describes the infectious compartment , shown that a solution to this formulation exists and that the solution is unique. We have then gone on to demonstrate an algorithm to approximate the solution to the nonlinear observational model. While it is standard these days for more complex and descriptive compartmental models to be analysed and used in practice, this study provides a first illustration of the analysis for the novel boundary value approach, which we intend to build upon to combat the issues outlined in §§2 and 3. It is not hard to reason that using two data points is not enough to provide accurate forecasting capabilities for the short and medium term, but this result combined with parameter estimation with more data points and a model selection approach will shed light on how one can use this approach.

Now that we have guaranteed that the solution to the observational model is unique, given some parameters, we next look to answer the questions of the identifiability of the parameters from the data and how to appropriately deal with noise in the data. That is to say, how much data do we need to be able to uniquely identify the parameters and as well as ? These are questions of particular importance as understanding what we can identify, and the accuracy of the estimations, helps us to understand the scenario-based forecasting capabilities of the mathematical model. Understanding this and therefore the underlying assumptions needed to produce forecasting results, given the limelight of compartmental models in recent history, will allow for a confident and trustworthy exchange of knowledge between mathematical modellers and those with research questions on infectious diseases, such as local government or healthcare providers.