A dynamic model for the COVID-19 with direct and indirect transmission pathways.

Two common transmission pathways for the spread of COVID-19 virus are direct and indirect. The direct pathway refers to the person-to-person transmission between susceptibles and infectious individuals. Infected individuals shed virus on the objects, and new infections arise through touching a contaminated surface; this refers to the indirect transmission pathway. We model the direct and indirect transmission pathways with a S A D O I R ode model. Our proposal explicitly includes compartments for the contaminated objects, susceptible individuals, asymptomatic infectious, detected infectious, and recovered individuals. We compute the basic reproduction number and epidemic growth rate of the model and determine how these fundamental quantities relate to the transmission rate of the pathways. We further study the relationship between the rate of loss of immunity and the occurrence of backward bifurcation. An efficient statistical framework is introduced to estimate the parameters of the model. We show the performance of the model in the simulation scenarios and the real data from the COVID-19 daily cases in South Korea.

INTRODUCTION

Coronavirus disease 2019 (COVID‐19) has now become an unprecedented challenge around the world, a disease that began in December 2019 in Wuhan, China, and quickly spread to all parts of the world. The World Health Organization declared the disease to be a global pandemic in March 2020. In addition to the high number of cases and deaths in the world, this pandemic has had dramatic and unprecedented effects on many aspects of the human lives such as economic and industrial activities, education, health, religion, recreation, and entertainment to name some. Since the first days of the disease, many efforts have been dedicated to helping to understand and analyze the dynamic nature of the spread of the virus and the influential factors. , , , , , ,

Coronavirus transmits through two transmission pathways. The direct pathway that is the infection occurs by droplet inhalation, i.e., aerosol pathway. Furthermore, infected individuals contaminate objects by touching, droplet shedding through coughing, sneezing, or exhaling and cause the pathogen to transfer from contaminated objects to hands and then to the mouth, nose, and eyes of individuals which causes new infectious cases; this refers to the indirect pathway. Some studies highlight the effect of the indirect transmission pathway in influenza and other infectious diseases; see previous works. , , , , ,

In this manuscript, we propose a multi‐compartment vector‐borne like model that incorporates both direct and indirect transmission pathways and the contaminated objects.

This paper is organized as follows. In Section 2, we present the model and prove the positivity and boundedness of the solutions. In Section 3, we compute the basic reproduction number and the initial epidemic growth rate of the model and study the effect of variation in the transmission rates of the two transmission pathways. In most of the models proposed in the literature for the COVID‐19,  , , , , , , transmission from the indirect route is not considered. Our study shows that with a specified mean asymptomatic infectious period, ignoring the indirect transmission pathway could lead to underestimating the basic reproduction number and overestimating the epidemic growth rate. In Section 4, we prove that if the rate of loss of immunity exceeds a certain value, backward bifurcation occurs, that leads to bistability and makes it more difficult to control the disease. Finally, in Sections 5, 6, and 7, we show the performance of the model in the simulation studies and apply the methodology to the COVID‐19 daily cases in South Korea.

THE MODEL FORMULATION

Our model has the following compartments: the class S for susceptible individuals; the class A that consists of two groups of individuals in the community, asymptomatic infectious individuals, i.e., infected individuals who have no symptoms of the disease and individuals with symptoms of the disease that are not detected by the healthcare system; the class D for infectious individuals who are detected by the healthcare system, i.e., confirmed cases; and the class of recovered individuals R. We assume that all detected cases are hospitalized until recovery or death. We also consider two compartments for the objects, clean objects and infectious, i.e., contaminated objects, .

Individuals in the class A spread the disease between the community through two routes: (I) direct route with the flow where is the probability of virus transmission through direct contact; and (II) indirect route, i.e., the transmission of infection by touching contaminated objects.

Objects are defined as surface and items in the environment that can be touched by an individual's hands. We exclude those objects which are normally never touched by humans, such as under the tables, car floor, and ceilings. These objects may not be homogeneously distributed around, for example, some of these objects may be more accessible, such as an elevator. Or some objects may be more likely to be contaminated, for example, objects in hospitals and other medical facilities. On the other hand, the type of surface and the temperature or the humidity of the surroundings could affect the effect of the objects on virus transmission. We ignore the heterogeneity in objects and their pollution and assume that individuals touch objects homogeneously and also pathogen agents are homogeneously distributed on objects.

Individuals in the compartment A transmit contamination to the clean objects by touching and shedding pathogen agent to the surface of the objects through coughing, sneezing, and exhaling at rate . We use the term for the indirect route flow from S to A where is the touching rate, with K the mean number of objects an individual touches in a unit of time, and the total number of objects in the environment. Furthermore, shows the fraction of contaminated objects, and is the probability of getting infected due to touching contaminated objects.

The fractions and are detected by the healthcare system, due to respective exacerbation of symptoms of the disease and identification measures, such as measuring fever in urban and rural passages. The fractions and represent the rate of recovery from the disease and the fraction of the recovered individuals lose their immunity and become susceptible again, is called the rate of loss of immunity.

The period of the disease may be long; hence, for the study of long term dynamics of the disease, we consider the recruitment Λ and the natural death rate in the model. Furthermore, we accommodate the disease‐induced mortality rate of m into the model. Finally, we consider as the removal rate of the pathogen agent from the objects or equivalently as the mean lifetime of the virus on the objects.

Based on the flow diagram of the model depicted in Figures 1 and 2, we define the following system of ordinary differential equations: with . Furthermore, we have the following system of equations for objects:

FIGURE 1

The flowchart of the individuals [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 2

The flowchart of the objects [Colour figure can be viewed at wileyonlinelibrary.com]

Let be constant; hence, we obtain the following equation:

In Lemma 2.1, we prove that the solutions are nonnegative.

If the initial conditions are nonnegative, i.e., , and , then all components of the solution ( in the system are nonnegative for all .

All components of the solution ( of the system are continuously differentiable. Furthermore, if all compartments have nonnegative initial conditions and if a compartment has zero value at time , then the its derivative is nonnegative. For example, if , , , , and , we get which implies ; hence, is nonnegative for all time . Next, assuming that , , , , and , we have which implies ; hence, is nonnegative for all time . The same discussion is true for other compartments and as mentioned in Ng and Gui, it implies the nonnegativity of all compartments at all time .

Boundedness is one of the basic properties of the solutions that we prove in Lemma 2.2.

For any nonnegative initial values, the total population and are bounded from above.

We have , and integration yields, and for all . This proves the boundedness of the solutions of system.

BASIC REPRODUCTION NUMBER AND EPIDEMIC GROWTH RATE

The system has the following DFE, i.e., disease free equilibrium point, . For the computation of the basic reproduction number, we use the linearization theorem.

The disease‐free steady state is locally asymptotically stable if and unstable if with

. The Jacobian matrix of the system at the point has the following form: that has the eigenvalues and the eigenvalues of the following submatrix: From the Routh‐Hurwitz criterion, it is well known that in a 2 × 2 matrix B, all eigenvalues have negative real parts if and only if and . Furthermore,

Now the relation can be rewritten as which implies the relation .

We define the direct and the indirect basic reproduction numbers by

Note that when there is no indirect transmission, i.e., , basic reproduction number reduces to . And the system (1) reduces to the following system:

In the early stages of the epidemic, growth is exponential and the initial epidemic growth rate is a measure of the severity of the epidemic, which is closely related to the basic reproduction number of the disease; see Chowell et al. and Ma. Here, we study the relation of these quantities in our model and the effect of variation in the direct and indirect transmission rates.

Given , a disease outbreak occurs and the Jacobian matrix at the DFE has a positive eigenvalue. Let denote the dominant eigenvalue of the Jacobian of the system at the disease‐free steady state , which is referred to as the initial outbreak growth rate; see Pinheiro et al. Calculating this eigenvalue yields the following formula: where with the additional property . Note that when there is not any indirect transmission i.e., , system (1) collapses to model (8), and .

In Wallinga and Lipsitch, it is studied how the combination of empirically observed growth rates and knowledge of the distribution of disease generation time estimates . In our model, the generation time consists of the mean asymptomatic infectious period and the virus lifetime in the object compartment. Let denote the mean time spent in asymptomatic infectious stage, and denote the mean virus lifetime on the surface of objects. After some algebra, (9) can be rewritten as the following relation between and :

When there is no indirect transmission, the relation (11) becomes the following relation in model: where . Estimation of basic reproduction number and initial epidemic growth rate from time series data is possible. Suppose that observed growth rate and observed basic reproduction number, i.e., and , are available from an incidence time series data. We compare the basic reproduction numbers and epidemic growth rates of the and models with the same underlying data set, in order to illustrate how the object compartment affects the estimation of and .

Suppose is assumed. Consider two set of parameters for the and systems, for which the epidemic growth rate in both models are the same as . Furthermore, assume the mean time spent in asymptomatic infectious stage, i.e., , be equal in both models and , denote the corresponding basic reproduction number, respectively. Then .

The relations (11) and (12) transform to and , respectively. By subtracting these two relations, we have

Now since , we obtain .

From a practical point of view, the above proposition states that when the mean asymptomatic infectious period is known, ignoring the indirect transmission route and fitting an model to an observed epidemic growth rate causes an underestimate of the basic reproduction number.

Let . Consider two sets of parameters for the and systems. Assuming the basic reproduction number in two models are the same as . Furthermore, assume the mean time spent in asymptomatic infectious stage, i.e., , be equal in both models and , denote the corresponding initial epidemic growth rates, respectively. Then .

The relations (11) and (12) transform to and , respectively. By subtracting these two relations, we have

Because , we obtain .

The above proposition implies that when the mean asymptomatic infectious period is known, deleting the indirect transmission route and fitting an model to an observed basic reproduction number lead to overestimating the epidemic growth rate of the model.

In the next proposition, we study the effect of variation in the contribution of direct and indirect transmission pathways on the disease dynamics, under fixed , by computing the total derivative of the epidemic growth rate with respect to .

Assume has a fixed value for a set of parameters in the system. And be the total derivative of the epidemic growth rate with respect to , respectively. Then and .

The quantity depends on and . Because we assume the basic reproduction number is fixed, and has the following relation: From relation (10) and the above relation, we get . Taking derivative from both sides of (11) with respect to , that yields

Because and , this shows . From similar computations, we get that shows .

The above proposition shows that, for fixed , increasing the contribution of the direct transmission results in a faster epidemic growth rate and in contrast increasing the contribution of the indirect transmission pathway results in a slower epidemic growth rate. As it has been mentioned in Tien and Earn for water‐borne diseases, disease transmission through the objects can be thought of as a delayed transmission route, as pathogens must first pass through the objects compartment before reaching a susceptible host. In fact, an infected individual has an effective infectious period, which takes into account both person‐person and person‐object‐person transmission. This effective infectious period increases as the relative contribution of the objects transmission increases. In order for to be fixed, the epidemic growth rate must decrease to compensate. This is an intuitive understanding of the above result.

Finally in this section, we study the effect of variation in the contribution of direct and indirect transmission pathways on the disease dynamics, under fixed epidemic growth rate .

Assume has a fixed value for a set of parameters in the system. And be the total derivative of the basic reproduction number with respect to , respectively. Then and .

The quantity depends on both and . Since we assume the epidemic growth rate is fixed, and are dependent together with and by differentiating from (9) with respect to , we have the following relation: Now by using (11) and since is constant, we have

The above proposition shows that, for fixed , increasing the contribution of the direct transmission reduces the basic reproduction number, and in contrast, increasing the contribution of the indirect transmission pathway increases the basic reproduction number.

BACKWARD BIFURCATION AND ENDEMIC EQUILIBRIUM POINTS

Generally in epidemiologic models, the spread of infection can be controlled by reducing to the region , if the initial size of all compartments of the model are in the basin of attraction of the DFE . But also, in some models, in the range , endemic equilibrium points may also exist, which shows that reducing to is not enough for eliminating the disease. In such models, it is said that backward bifurcation occurs; see Martcheva. In this section, we study the occurrence of backward bifurcation in our model.

A standard method for studying the occurrence of backward bifurcation is to use the Castillo‐Chavez and Song theorem, that was obtained by the center manifold theorem. The following result is obtained by directly applying theorem 4.1 in Castillo‐Chavez and Song.

Backward bifurcation occurs in the model at , when where .

First, we convert the system variables as follows: , , , , and , so that our system transforms into the following system:

We consider as the bifurcation parameter. Hence, when , we have , and if and only if and if and only if . Now if , we have which implies that 0 is a simple eigenvalue of , and other eigenvalues are negative. Let be the right eigenvector of associated with the eigenvalue 0, i.e., , computation yields: Similarly, we obtain the following vector , as the left eigenvector of associated with the eigenvalue 0, i.e., : Now as it is proved in theorem 4.1. of Castillo‐Chavez and Song, if the bifurcation quantities and are both positive, then backward bifurcation occurs in the system, and we have Furthermore, implies

The above theorem states that when , backward bifurcation cannot occur. Hence, the disease‐free equilibrium point can be global asymptotic stable, which ensures that the elimination of the disease is independent of the initial size of the compartments of the model.

When , the disease‐free equilibrium point is globally asymptotically stable for .

We use the criterion introduced in Castillo‐Chavez et al. for global stability. We rewrite our system in the form (3.1) in Castillo‐Chavez et al. as the following form: with and Now it is clear that and is global asymptotic stable in ; hence, is global asymptotic stable.

Finally, we compute the endemic equilibrium points of the model in term of parameters. Let us denote the endemic equilibrium point by . From (1), (3), and , we have where and is the positive root of the following quadratic equation: with the following coefficients:

PARAMETER ESTIMATION

Estimation of the model parameters and initial population size is one of the most important parts of studying epidemic mathematical models. The model introduced in (1) is dependent on both unknown parameter values and the initial population size. Usually in such models, the model parameters are estimated by assuming the initial population size and the estimation of the model with data from the observed COVID‐19 cases. In this section, we are going to introduce a statistical framework to estimate the important parameters of the model in (1) and then evaluate the estimated model with real data set. Our main motivation for using this parameter estimation approach is to have a partially observed system feature. We will apply this critical feature of our approach in daily observed data COVID‐19 for the South Korea data set. Below, we introduce our parameter estimation model.

Generalized Tikhonov regularization for ODE estimation

In the previous sections, we proposed an ODE system for the spread of COVID‐19 with both direct and indirect transmission. Our system contains several parameters, such as general mortality rate, the mortality rate due to COVID‐19, and the recovery rate of confirmed cases, can be retrieved from medical references, and general information. On the other hand, some parameters of the system such as direct and indirect transmission rate, object contamination rate, exacerbation rate, and detection rate can be estimated by fitting the model with the daily data of the detected, i.e., confirmed cases. However, a problem in the fitting of the proposed ODE system is the partial observation of the system; i.e., we have data only in the state D of the system. ,

Recently, Vujacic et al. show that the generalized Tikhonov regularization (GTR) is a general structure for parameter estimation in ordinary differential equations. They showed that generalized Tikhonov regularization can be applied to partial measured systems. Partial measurement systems refer to a system in which common states or specific variables are expensive in the model that cannot be measured or when it is not possible to measure or view data for a particular state for various reasons. For example, in our system in (1), there is only observational data from state D, and for other states in the model, it is not possible to collect observational data. This led us to use the generalized Tikhonov regularization as a tool for parameter estimation in this section.

Compared to other approaches, the parameters estimated from this method have minimal bias and variance and are as good as other approaches, and compared to other methods such as nonlinear least square (NLS), which adds solution approximation noise to data noise, it performs better. The generalized Tikhonov function for the general form of ODE is as follows: where the functions , , and are defined in Vujacic et al. as Objective function, Stabilizing function, and Similarity function, respectively.

This approach incorporates a collocation method and deviations from differential equation models for estimation.

We have written the vector form for showing the more general notation as and refer to the components of by for . In this short‐hand, we rewrite Equations (1) and (3) as where vector shows all the unknown parameters in model as .

Our aim is to infer given n noisy observations of the solution of the differential equation (1), , where where are the d‐dimensional column vectors (here ) of measurement errors at time and .

The regularized solution is found by optimizing (25) over function space parametrized by . We take two steps to estimate the parameters. As a first step, we assume parameter is fixed, then any class of was estimated by as the equal dimensional function space with basis functions { belong to dimension . With using the access of , , by using data from (26), we have where and . We have used B‐splines as common basis functions to achieve a sequence of spaces whose union is dense in (.

In the latter step, we optimize (25) with respect to over : As mentioned above, Vujacic et al. introduce this approach as a general method to estimate parameters in ODE models which can find more details there.

NUMERICAL STUDIES

The performance of the proposed model in Equation (1) on simulated and real data is studied in this section. To reduce the complexity of the parameter estimation for the saturated model with 12 unknown parameters (introduced in Table 1), we focus on a set of five parameters selected by partial rank correlation coefficient (PRCC), namely, the direct transmission rate , indirect transmission rate , rate of objects contamination , rate of exacerbation of symptoms of the disease , and identification rate of infectious individuals . The remaining seven parameters as well as the initial conditions are chosen to be aligned with the literature , , and shown in Table 2.

TABLE 1

Description of variables and parameters of the model. The appropriate location for the table 1 is in Sect 2

Symbol	Description
S	Susceptible individuals
A	Infectious individuals without symptoms + individuals with symptoms that are not detected by the healthcare system
D	Infectious individuals who are detected by the healthcare system, i.e., confirmed cases
R	Recovered individuals
O I	Contaminated objects
Λ	Recruitment rate into the susceptible population
μ	Natural death rate
m	Death rate due to COVID‐19
β 1	Direct transmission rate
β 3	Indirect transmission rate
δ	Object contamination rate
δ 1	Rate of exacerbation of symptoms of the disease
γ	Loss of immunity rate
γ 1	Rate of recovery of confirmed cases
ν	Identification rate of infectious individuals
ξ	Rate of recovery of asymptomatic infectious individuals
1ξ_1	Mean lifetime of virus on the objects
K	Mean number of objects an individual touches in a unit of time
N O	Total number of objects in the environment

TABLE 2

The initial values and the conditions for the seven parameters that are determined to be less sensitive by PRCC

Parameter	μ	m	ξ	γ 1	Λ	ξ 1	γ	N 0
Value	10−6	0.02	0.2	0.2	50	0.1428	10−10	109
Initial condition	S(0)	A(0)	D(0)	R(0)	O I (0)
Value	53 000 000	190	19	0	20 000

Note: The values are set to be align with the literature.

Simulation study

We show the performance of the generalized Tikhonov regularization (GTR) to estimate the selected five parameters, (, in simulation studies. The method is applied to the different scenarios with different sample sizes n and data generated as in Equation (26). Each simulation is repeated 100 times, and the results showed are empirical means across those runs. We first generate the solutions of over the time interval based on initial condition in Table 2 and , , , , . We add Gaussian white noise ) to each solution. Sample sizes are varied according to . The initial guess in the optimization of the is set to the true value of the parameters. For the regularization term or penalty parameter, we assume exponential penalties and . We utilized the third‐order B‐splines and nonlinear optimizer nlminb in R for the inner and outer optimizations.

The results are presented in Table 3. As GTR suggests, we see that as the sample size increases as well as the regulatory parameter decreases, the convergence of the estimated values to the true values improves.

TABLE 3

The table of true and estimated values from GTR for the simulation studies

Penalty parameter	n	β 1	β 3	δ	δδ 1	ν
Gaussian error, σ=0.1
η=10^−2	10	2.55 × 10−8	6.3 × 10−15	1.93 × 10−8	5.266	3.115
	100	1.89 × 10−7	2.745 × 10−10	4.53 × 10−7	2.021	1.724
	10 000	2.79 × 10−7	4.3 × 10−11	2.93 × 10−6	0.866	0.515
η=10^−3	10	1.22 × 10−7	5.7 × 10−14	5.44 × 10−7	1.09	0.527
	100	1.94 × 10−7	3.08 × 10−12	4.01 × 10−6	0.09	0.189
	10 000	2.95 × 10−7	3.83 × 10−11	3.07 × 10−6	0.166	0.116
Gaussian error, σ=0.01
η=10^−2	10	2.67 × 10−8	5.1 × 10−15	1.07 × 10−8	4.266	2.115
	100	1.94 × 10−7	1.32 × 10−12	3.12 × 10−7	1.812	1.061
	10 000	2.98 × 10−7	3.9 × 10−11	2.97 × 10−6	0.139	0.112
η=10^−3	10	1.04 × 10−7	2.06 × 10−12	1.04 × 10−7	0.217	0.223
	100	2.24 × 10−7	2.98 × 10−11	3.89 × 10−6	0.101	0.118
	10 000	3.09 × 10−7	3.11 × 10−11	3.12 × 10−6	0.131	0.124
True parameter value		3.10 × 10−7	3.12 × 10−11	3.105 ×  10−6	0.132	0.125

Note: This table shows the performance of the model for selected five parameters under two scenarios with the sample size as well as penalty value .

The same simulation scenario is repeated with higher noise level, , and the results are shown in complementary Figures 3 and 4. The simulated data are represented by the circles. The curves are obtained by solving the model (1) for true parameters (red line) and the curves are obtained by solving the model (1) for the estimated parameter (blue line). These figures show that our method converges to the true curve (red line in the figures) for D, R, and O states, however, with a significantly high variation for the rest of the parameters.

FIGURE 3

The performance of the GTR to estimate the selected five parameters for the simulated data with . The circles show the generated data, the red line represents the solution of the model (1) based on true parameter values with , and the blue line represents the solution of the model (1) based on estimated parameter values [Colour figure can be viewed at wileyonlinelibrary.com]

FIGURE 4

The performance of the GTR to estimate the selected five parameters for the simulated data with . The circles show the generated data, the red line represents the solution of the model (1) based on true parameter values with , and the blue line represents the solution of the model (1) based on estimated parameter values [Colour figure can be viewed at wileyonlinelibrary.com]

Real data application

In this section, we show the application of the GTR to the daily COVID‐19 cases in South Korea collected from https://www.worldometers.info/coronavirus/. These data are examples of a partially observed system. In other words, the states are not fully known (here, only state D is known) to the model, and as a result, we estimate the solution of states and the parameters with data only from state D.

The outcome of the parameter estimation is reported in Table 4 for different penalty parameters. Figure 5 shows the fitted curve to the data. The observations are displayed with circles while the solution based on the GTR estimate for four states solution of model (1) is plotted with a different color line.

TABLE 4

The fitting of GTR to the COVID‐19 daily infections in South Korea

Parameter	β 1	β 3	δ	δ 1	ν
Estimated value	1.1905 × 10−5	2.9056 × 10−9	8.9221 × 10−6	0.1186	0.3259

FIGURE 5

State D in the model (1), the circles represent the South Korean data, the blue line represents the estimated solution. Other colors represent estimated solution for the other states [Colour figure can be viewed at wileyonlinelibrary.com]

The results of estimating the model parameters are reported in Table 4. As we have already mentioned, these results are obtained only by observing daily data from state D, and all the sensitive parameters in the model are estimated.

CONCLUSION

A model is developed for the spread of COVID‐19 incorporating both direct and indirect transmission pathways. The basic reproduction number and the epidemic growth rate of the model are obtained and the effect of variation in the contribution of various routes of transmission on these quantities is studied. Our study shows that with a known mean asymptomatic infectious period, neglecting the indirect transmission route and fitting the model to the observed epidemic growth rate leads to overestimating the basic reproduction number of the model. And in contrast, we show that neglecting the indirect transmission route and fitting model to the observed basic reproduction number leads to overestimating the epidemic growth rate in the model. Furthermore, it is proven that if the rate of loss of immunity exceeds a certain value, backward bifurcation occurs, that leads to bistability and makes it more difficult to control the disease.

We further have used a generalized Tikhonov regularization method to estimate the unknown parameters from observational data. In simulation studies, we have shown that the dynamic model can be identifiable from noisy data when the data are generated for all five states (class). The simulation results confirmed our theoretical identifiability analysis. Although our simulation study has focused on the 5D dynamic model, in reality, we normally encounter partially observed data. Therefore, the 5D dynamic model was fitted to daily data of COVID‐19 for South Korea detected cases only for class D consist of detected infectious individuals i.e., confirmed cases of COVID‐19. Our proposed approaches work for this case as well and the basic reproduction number and the epidemic growth rate were estimated.