Stability and bifurcation analysis of a nested multi-scale model for COVID-19 viral infection

Bishal Chhetri; Krishna Kiran Vamsi Dasu

doi:10.1515/cmb-2024-0006

Article Open Access

Stability and bifurcation analysis of a nested multi-scale model for COVID-19 viral infection

Bishal Chhetri and Krishna Kiran Vamsi Dasu

Published/Copyright: June 21, 2024

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Computational and Mathematical Biophysics Volume 12 Issue 1

Abstract

In this study, a nested multi-scale mathematical model is proposed and studied for COVID-19 viral infection. The well-posedness of the multi-scale model is discussed, followed by the stability analysis of the equilibrium points. The proposed model admits two equilibrium solutions: the infection-free equilibrium and the infected equilibrium. The infection-free equilibrium point is shown to be globally asymptotically stable provided that the value of the basic reproduction number is less than unity. When the value of the basic reproduction number exceeds unity, a unique infected equilibrium exists, and the system is found to undergo a forward (trans-critical) bifurcation. Two-parameter heat plots are performed to identify the parameter combinations for the stability of the admitted equilibrium points. From the sensitivity analysis, the transmission rate, infection rate, and birth rate were found to be most sensitive to the basic reproduction number. The influence of the within-host sub-model parameters on the between-host sub-model variables is numerically illustrated. The spread of infection in a community is shown to be influenced by the within-host sub-model parameters, such as the production of viral particles, the clearance rate of infected cells by the immune system, and the clearance rate of viral particles by the immune system. The comparative effectiveness of three health interventions (antiviral drugs, immunomodulators, and generalized social distancing) for COVID-19 infection is examined using the reproductive number as an indicator of the effectiveness of these interventions. The results from the comparative effectiveness study suggest that a combined strategy of antiviral drugs, immunomodulators, and generalized social distancing would be the best strategy to implement to contain the spread of infection in the community.

Keywords: multi-scale model; well-posedness; elasticity; comparative effectiveness; heat plot

MSC 2010: 37; 58; 65

1 Introduction

Infectious diseases with their host-parasite interactions represent a complex system involving components and processes that operate at different levels of time, space, and biological organization [19]. Infectious disease dynamics involve thousands of molecules, complicated interactions of the pathogen with the host and also with the environment. Therefore, to understand disease dynamics in-depth, one must understand the processes at these different levels of manifestation of infectious diseases, and the interactions between these levels. The two most studied scales in the transmission of infectious disease systems are the within-host and between-host scales. Although these models have provided useful information about virus-cell interactions at the level of a single host and disease transmission between hosts, disease modelers see the need to integrate or link these models to obtain a better understanding of the dynamics of infectious disease spread.

Several studies on infectious diseases have shown that disease transmission depends on within-host processes [2,14,21,23,24,31,37]. In particular, these studies have shown that the transmission potential of an infected host depends on the pathogen load and increases with increasing pathogen load. The different functional relationships between host infectivity and pathogen load have been identified. Studies of dengue and human immunodeficiency virus (HIV) transmission revealed a sigmoidal functional relationship between host infectivity and the pathogen load at the within-host level [31,37]. Studies that discuss similar functional relationships between host infectivity and pathogen load at the within-host scale include [23] for malaria and [24] for human T lymphotropic virus, [32] for visceral leishmaniasis and chikungunya, and [46] for hepatitis B. Over the past decade, the importance of multi-scale modeling of infectious disease dynamics has been increasingly recognized, and many mathematical models linking within-host and between-host processes have been developed. The multi-scale models for influenza infections can be found in [18,22]. In [17], the within-host model is nested within a broader epidemiological model by linking the transmission rate of the infection or the additional host mortality rate to the dynamics of the within-host model. In studies such as [2,14], a nested multi-scale model (NMSM) is developed that considers the transmission rate at the intermediate host level as a linear function of viral load. A multi-scale model is developed and used to evaluate the effectiveness of health interventions acting at different time scales in [16]. Using a similar dynamic approach, there is considerable interest in linking disease dynamics at within-host and between-host levels, such as in the study of HIV and hepatitis C virus [11,30]. The multi-scale models for influenza infections can be found in [18,22]. Multi-scale mathematical modeling studies of infectious diseases have expanded the landscape of disease modeling and have the potential to more accurately describe disease dynamics [33]. A thorough understanding of the infectious disease transmission requires knowledge of the processes at different levels of the infectious disease and of the interplay between these levels. To obtain a clear idea of the dynamics of the disease, different time scale models need to be integrated [2,14,21].

Multi-scale models of the infectious disease systems integrate the within-host scale and the between-host scale. There are five main different categories of multi-scale models that can be developed at different levels of organization of an infectious disease system, which are: individual-based multi-scale models (IMSMs), NMSMs, embedded multi-scale models (EMSMs), hybrid multi-scale models (HMSMs), and coupled multi-scale models (CMSMs) [15]. In IMSMs, the within-host sub-model is used to describe the entire infectious disease system across both the within-host scale and between-host scale. In case of NMSMs, there is only unidirectional flow of information (only from within-host scale sub-model to between-host scale sub-model). There is never, any reciprocal feedback from the between-host scale sub-model back down to the within-host scale sub-model. In contrast, EMSMs is both a top-down and bottom-up modeling approach where both the within-host scale sub-model and the between-host sub-model influence each other. In all the aforementioned three categories of multi-scale models, the within-host scale and the between-host scale belonged to the same domain and is modeled in a homogeneous way using the same formalism. But in case of HMSMs, the within-host scale and the between-host scale may no longer be modeled in a homogeneous way. The HMSMs are formed based on the multi-domain integration framework. The details of these categories can be found in [15].

Multi-scale modeling of COVID-19 disease is still in its infancy, and efforts are being made in this direction. Some of the available articles on multi-scale modeling can be found in [1,4,36,45]. A multi-scale model would be extremely helpful in understanding the spread of COVID-19 infection and evaluating the efficacy of the interventions not only at the individual level but also at the population level. Therefore, in this study, we develop an NMSM that integrates within-host and between-host sub-models. As the importance of multi-scale modeling in disease dynamics is increasingly recognized, we believe that our study contributes to the growing knowledge on multi-scale modeling of COVID-19 disease.

This article is organized as follows: in Section 2, we develop an NMSM, study the stability and bifurcation analysis, and numerically illustrate the theoretical results obtained. We also analyze the sensitivity of R 0 with respect to the model parameters and present two-parameter heat plots to identify the regions of stability. In Section 3, we conduct a comparative effectiveness study to evaluate the effectiveness of public health interventions for COVID-19. In Section 4, we provide the discussion and conclusions.

2 Multi-scale model

The multi-scale model that we develop describes COVID-19 disease dynamics across two different time scales, i.e., the within-host scale and the between-host scale. The within-host scale model tracks the density of the pathogens and also the state of the host’s defense mechanisms, while the between-host scale model is concerned with the transmission of the disease within the host population. The multi-scale model consists of seven variables, namely, susceptible cells ( U ) , infected cells ( U * ) , viral load ( V ) , and four between-host variables, namely, susceptible human population ( S ) , infected human population ( I ) , exposed human population ( E ) , and recovered human population ( R ) . We work with the following assumptions for the multi-scale model as in [16,36].

The dynamics of the within-host scale variables are assumed to occur at a fast time scale s so that U = U ( s ) , U * = U * ( s ) , and V = V ( s ) , while the dynamics of the between-host scale variables is assumed to occur at a slow time scale t so that S = S ( t ) , I = I ( t ) , E = E ( t ) , and R = R ( t ) . Because the infection obtains cleared within a few days, the value of V ( s ) is non-zero only for a short period of time. Due to this, the within-host dynamics is considered only for a short period of time. But because the infection remains for longer time at the population level, the between-host dynamics is considered for a longer period of time.
The transmission rate β and disease-induced death rate d at the between-host scale are assumed to be functions of viral load.
Immune response is captured in the model through the parameters x and y ; these parameters, respectively, denote the combined rates at which the infected cells and the virus are cleared by the release of cytokines and chemokines such as IL-6 TNF- α , INF- α , CCL5, CXCL8 and CXCL10 [10].

Based on the aforementioned assumptions, the multi-scale model is described by the following system of differential equations:

(2.1) d U d s = ω − k U ( s ) V ( s ) − μ c U ( s ) ,

(2.2) d U * d s = k U ( s ) V ( s ) − x U * ( s ) − μ c U * ( s ) ,

(2.3) d V d s = α U * ( s ) − y V ( s ) − μ v V ( s ) ,

(2.4) d S d t = Λ − β ( V ( s ) ) S ( t ) I ( t ) − μ S ( t ) ,

(2.5) d E d t = β ( V ( s ) ) S ( t ) I ( t ) − ( μ + π + γ 1 ) E ( t ) ,

(2.6) d I d t = π E ( t ) − ( μ + γ 2 ) I ( t ) − d ( V ( s ) ) I ( t ) ,

(2.7) d R d t = γ 1 E ( t ) + γ 2 I ( t ) − μ R ( t ) .

The within-host model (2.1)–(2.3) considered in this study is studied in detail by the authors in [10]. The meaning of each of the variables and parameters of the model (2.1)–(2.7) is given in Table 1.

Table 1

Meanings of the parameters

Parameters	Biological meaning
ω	Natural birth rate of cells
Λ	Birth rate of human population
α	Burst rate of the virus
μ	Natural death rate of human population
μ c	Natural death rate of cells
μ v	Natural death rate of virus
π	Infection rate of exposed population
k	Infection rate of susceptible cell
γ 1 , γ 2	Recovery rate of the exposed and infected human population
β	Transmission rate of the susceptible human population
x	Rate at which the infected cells are cleared by the individual immune response
y	Rate at which the virus particles are cleared by the individual immune response
d	Death rate of the infected human population

Based on the experimental observations of the impact of viral load on disease transmission [25] and disease-induced deaths, the linking of within-host and between-host sub-models is modeled by means of β = β ( V ( s ) ) and d = d ( V ( s ) ) . Though the exact functional relationship between viral load, transmission rate, and disease-induced death rate is not known [2], as in [2,14], we considered a linear form for the coupling functions β and d : β ( V ( s ) ) = β V ( s ) and d ( V ( s ) ) = d V ( s ) . We also observe that the first three equations of the between-host SEIR sub-model are independent of the recovered population R ( t ) . In similar lines to the studies done in [26,29], without loss of generality, we omit the last equation of the between-host SEIR sub-model. Taking β and d as linear functions of viral load and omitting the equation for R ( t ) in the between-host SEIR sub-model, the final multi-scale model for COVID-19 disease dynamics is given by the following set of equations:

(2.8) d U d s = ω − k U ( s ) V ( s ) − μ c U ( s ) ,

(2.9) d U * d s = k U ( s ) V ( s ) − x U * ( s ) − μ c U * ( s ) ,

(2.10) d V d s = α U * ( s ) − y V ( s ) − μ v V ( s ) ,

(2.11) d S d t = Λ − β V ( s ) S ( t ) I ( t ) − μ S ( t ) ,

(2.12) d E d t = β V ( s ) S ( t ) I ( t ) − ( μ + π + γ 1 ) E ( t ) ,

(2.13) d I d t = π E ( t ) − ( μ + γ 2 ) I ( t ) − d V ( s ) I ( t ) .

2.1 Reduced multi-scale model for COVID-19 dynamics

The multi-scale model (2.8)–(2.13) is not that easy to analyze. The problem arises from the discrepancy between the time scales. The viral load V ( s ) at the within-host scale is a function of s , while the dynamics of the between-host variables are assumed to occur at a slow time scale t . To overcome these types of difficulties in multi-scale mathematical modeling studies of infectious diseases, a quantity called area under the viral load is defined using the within-host scale submodel [16,36]. In similar lines to the simplification done in [16,36], we simplify the multi-scale model (2.8)–(2.13) by defining the area under the viral load curve and use this quantity as a proxy for measuring the infectiousness of an individual host. We consider the within-host scale submodel from the multi-scale model (2.8)–(2.13) to derive the expression for the area under the viral load curve. From (2.8)–(2.13), the within-host submodel is given by the following system of differential equations:

(2.14) d U d s = ω − k U ( s ) V ( s ) − μ c U ( s ) ,

(2.15) d U * d s = k U ( s ) V ( s ) − x U * ( s ) − μ c U * ( s ) ,

(2.16) d V d s = α U * ( s ) − y V ( s ) − μ v V ( s ) .

The viral load V ( s ) at the within-host scale provides the link between the within-host scale and the between-host scale. We use the within-host scale sub-model given by the model system (2.14)–(2.16) to estimate the individual host infectiousness. To derive an expression for the area under the viral load curve, we use the within-host sub-model (2.14)–(2.16). We denote the area under the viral load curve by N h . The quantity N h measures the total amount of COVID-19 virus produced by an infected human throughout the period of host infectivity, and thus can be considered a proxy for individual host infectivity [20]. Let d 1 and d 2 be the times at which the viral load takes on the value of the detection limit at the beginning and end of infection. Then, d 2 − d 1 can be considered as the duration of host infectivity. The amount of virus produced N h , during the entire period of host infectiousness can be calculated from the within-host scale sub-model using the following formula:

N h = ∫ d 1 d 2 V ( s ) d s .

We obtain the expression for N h in similar lines to the method discussed in [16]. Integrating equation (2.16) from d 1 to d 2 and using the definition of N h , we obtain the following expressions:

V ( d 2 ) − V ( d 1 ) = α ∫ d 1 d 2 U * ( s ) d s − ( y + μ v ) ∫ d 1 d 2 V ( s ) d s .

Simplifying the aforementioned expression, we obtain

(2.17) N h = α ∫ d 1 d 2 U * d s + V ( d 1 ) − V ( d 2 ) y + μ v ,

where U * is the number of infected cells. The area under the viral load curve N h given by equation (2.17) can be numerically calculated, and it gives a fixed value for the chosen set of parameter values. Now, the within-host scale sub-model (2.14)–(2.16) is reduced to a composed parameter N h given by equation (2.17). This N h replaces V ( s ) in the between-host scale sub-model as follows:

(2.18) d S d t = Λ − β N h S ( t ) I ( t ) − μ S ( t ) ,

(2.19) d E d t = β N h S ( t ) I ( t ) − ( μ + π + γ 1 ) E ( t ) ,

(2.20) d I d t = π E ( t ) − ( μ + γ 2 ) I ( t ) − d N h I ( t ) .

We note here that the entire within-host scale submodel is reduced to a single parameter N h that depends on the within-host scale parameters. The reduced, simplified multi-scale model (2.18)–(2.20) is now on a single time scale t and is easier to analyze. The reduced multi-scale model (2.18)–(2.20) is categorized as an NMSM according to the categorization of multi-scale models for infectious diseases [15] and is unidirectional, i.e., the within-host submodel influences the between-host submodel without mutual feedback. We study the impact of within-host dynamics on between-host dynamics via the composite parameter N h .

2.2 Well-posedness of the model

The existence, the positivity, and the boundedness of the solutions of the proposed models (2.18)–(2.20) need to be proved to ensure that the model has a mathematical and biological meaning.

2.2.1 Positivity and boundedness

Positivity and boundedness of the solutions of Systems (2.18)–(2.20) are proved in similar lines to the method discussed in [3,29].

Positivity ̲ :

Lemma 2.1

Let S ( 0 ) ≥ 0 , E ( 0 ) ≥ 0 , and I ( 0 ) ≥ 0 , and then, the solutions S ( t ) , E ( t ) , and I ( t ) of Systems (2.18)–(2.20) are non-negative for all t ≥ 0 .

Proof

Let x ( t ) = ( S ( t ) , E ( t ) , and I ( t ) ) be the solution of the system with the initial condition x 0 = ( S ( 0 ) , E ( 0 ) , and I ( 0 ) ) . We assume that all the parameters of the model are positive. Because the quantity N h measures the total amount of COVID-19 virus produced by an infected human, we consider N h to be a positive parameter.

Now, using the continuity of solution, for all of S ( t ) , I ( t ) , and R ( t ) that have positive initial values at t = 0 , we have the existence of the solution on an interval ( 0 , t 0 ) such that S ( t ) ≥ 0 , I ( t ) ≥ 0 , and R ( t ) ≥ 0 for 0 < t < t 0 .

Now, let us say that we have t = t 1 ≥ 0 such that S ( t 1 ) = 0 and E ( t 1 ) ≥ 0 , I ( t 1 ) ≥ 0 . Then, from System (2.18), we find that at t = t 1

d S d t t = t 1 = Λ ≥ 0 .

This positive derivative ensures that at any time, the solution reaches the axis, its derivative increases, and the function S ( t ) does not cross to negative part. This shows that S ( t ) remains non-negative for all t ≥ 0 .

Similarly, using the same analysis, we can show that

d E d t t = t 1 = β N h S ( t 1 ) I ( t 1 ) ≥ 0 ,

d I d t t = t 1 = π E ( t 1 ) ≥ 0 .

Thus, we find that x ( t ) never crosses the axes S = 0 , I = 0 , and R = 0 when it touches them. Thus, for any positive initial conditions, the solutions S ( t ) , E ( t ) , and I ( t ) of Systems (2.18)–(2.20) remain non-negative for all t ≥ 0 .

Boundedness ̲ :

Let N ( t ) = S ( t ) + E ( t ) + I ( t )

Now,

d N d t = d S d t + d E d t + d I d t = Λ − μ ( S ( t ) + E ( t ) + I ( t ) ) − γ 1 E − γ 2 I − d N h ≤ Λ − μ N ( t ) .

Therefore,

d N d t + μ N ( t ) ≤ Λ .

The integrating factor is given by e μ t . Therefore, after integration, we obtain

N ( t ) ≤ Λ μ + c e − μ t , where c is a constant. Now, as t → ∞ , we obtain

lim sup N(t) ≤ Λ μ .

Thus, we have shown that the solutions of Systems (2.18)–(2.20) are ultimately bounded. Therefore, the biologically feasible region of Systems (2.18)–(2.20) is given by the following set:

□ Ω = ( S ( t ) , E ( t ) , I ( t ) ) ∈ R + 3 : S ( t ) + E ( t ) + I ( t ) ≤ Λ μ , t ≥ 0 .

We summarize the aforementioned discussion on boundedness by the following lemma.

Lemma 2.2

The set

Ω = ( S ( t ) , E ( t ) , I ( t ) ) ∈ R + 3 : S ( t ) + E ( t ) + I ( t ) ≤ Λ μ , t ≥ 0

is a positive invariant and an attracting set for Systems (2.18)–(2.20).

2.2.2 Existence and uniqueness of solution

For the general first-order ordinary differential equation of the form

(2.21) x ˙ = f ( t , x ) , x ( t 0 ) = x 0 .

One would have interest in knowing the answers to the following questions:

Under what conditions, solution exists for System (2.21)?
Under what conditions, unique solution exists for System (2.21)?

We use Theorem A.1 from the appendix to establish the existence and uniqueness of the solution for our susceptible, exposed and infected models (2.18)–(2.20).

Let D denote the domain:

∣ t − t 0 ∣ ≤ a , ∣ ∣ x − x 0 ∣ ∣ ≤ b , x = ( x 1 , x 2 , … , x n ) , x 0 = ( x 10 , … , x n 0 ) ,

and f ( t , x ) satisfies the Lipschitz condition:

(2.22) ∣ ∣ f ( t , x 2 ) − f ( t , x 1 ) ∣ ∣ ≤ k ∣ ∣ x 2 − x 1 ∣ ∣ .

Then, we have the following theorem.

Theorem 2.1

(Existence of solution) Let D be the domain defined above such that (2.22) holds. Then, there exists a unique solution of model system of equations (2.18)–(2.20), which is bounded in the domain D.

Proof

Let

(2.23) f 1 = Λ − β N h S ( t ) I ( t ) − μ S ( t ) ,

(2.24) f 2 = β N h S ( t ) I ( t ) − ( μ + π + γ 1 ) E ( t ) ,

(2.25) f 3 = π E ( t ) − ( μ + γ 2 ) I ( t ) − d N h I ( t ) .

We will show that

∂ f i ∂ x j , i , j = 1 , 2 , … , n ,

is continuous and bounded in the domain D.

From equation (2.23), we have

∂ f 1 ∂ S = − β N h I − μ , ∂ f 1 ∂ S = ∣ − β N h I − μ ∣ < ∞ , ∂ f 1 ∂ E = 0 , ∂ f 1 ∂ E < ∞ , ∂ f 1 ∂ I = − β N h S , ∂ f 1 ∂ I = ∣ − β N h S ∣ < ∞ .

Similarly, from equation (2.24), we have

∂ f 2 ∂ S = β N h I , ∂ f 2 ∂ S = ∣ β N h I ∣ < ∞ , ∂ f 2 ∂ E = − μ − π − γ 1 , ∂ f 2 ∂ E = ∣ − ( μ + π + γ 1 ) ∣ < ∞ , ∂ f 2 ∂ I = β N h S , ∂ f 2 ∂ I = ∣ β N h S ∣ < ∞ .

Finally, from (2.25), we have

∂ f 3 ∂ S = 0 , ∂ f 3 ∂ S < ∞ , ∂ f 3 ∂ E = π , ∂ f 3 ∂ E = ∣ π ∣ < ∞ , ∂ f 3 ∂ I = − ( μ + γ 2 + d N h ) , ∂ f 3 ∂ I = ∣ − ( μ + γ 2 + d N h ) ∣ < ∞ .

Hence, we have shown that all the partial derivatives are continuous and bounded. Therefore, Lipschitz condition (2.22) is satisfied. Hence, by Theorem 2.1, there exists a unique solution of Systems (2.18)–(2.20) in the region D .□

2.3 Stability analysis

The basic reproduction number denoted by R 0 that gives the average number of secondary cases per primary case is calculated using the next-generation matrix method [13], and the expression for R 0 for Systems (2.18)–(2.20) is given by

(2.26) R 0 = β N h π Λ μ ( μ + π + γ 1 ) ( μ + γ 2 + d N h )

Systems (2.18)–(2.20) admit two equilibria, namely, the infection-free equilibrium E 0 = Λ μ , 0 , 0 and the infected equilibrium E 1 = ( S * , E * , I * ) , where

S * = Λ ( β N h I * + μ ) ,

E * = ( μ + γ 2 + d N h ) I * π ,

I * = μ ( R 0 − 1 ) β N h .

Since negative population does not make sense, the existence condition for the infected equilibrium point E 1 is that R 0 > 1 .

2.3.1 Stability analysis of E 0

We analyze the stability of equilibrium points E 0 . This is done based on the nature of the eigenvalues of the Jacobian matrix evaluated at E 0 .

The Jacobian matrix of Systems (2.18)–(2.20) at the infection-free equilibrium E 0 is given by

J E 0 = − μ 0 − β N h Λ μ 0 − ( μ + π + γ 1 ) β N h Λ μ 0 π − ( μ + γ 2 + d N h ) .

The characteristic equation of J E 0 is given by

(2.27) ( − μ − λ ) ( λ 2 + ( 2 μ + π + γ 1 + γ 2 + d N h ) λ − ( R 0 − 1 ) ( μ + π + γ 1 ) ( μ + γ 2 + d N h ) ) .

One of the eigenvalues of characteristic equation ( 2.27 ) is − μ and the other two are the roots of the following quadratic equation:

(2.28) λ 2 + ( 2 μ + π + γ 1 + γ 2 + d N h ) λ − ( R 0 − 1 ) ( μ + π + γ 1 ) ( μ + γ 2 + d N h ) .

We see that when R 0 < 1 , both the eigenvalues of equation (2.28) are negative. Therefore, E 0 remains locally asymptotically stable whenever R 0 < 1 . When R 0 > 1 , the quadratic equation ( 2.28 ) has one positive and one negative root. Therefore, the characteristic equation (2.27) has two negative and one positive root. Hence, E 0 becomes unstable in this case. We summarize the local asymptotic stability of infection-free equilibrium E 0 in the following theorem.

Theorem 2.2

The infection-free equilibrium point E 0 of Systems (2.18)–(2.20) is locally asymptotically stable provided R 0 < 1 . If R 0 crosses unity, E 0 loses its stability and becomes unstable.

Global stability of E 0 ̲

To establish the global stability of the infection-free equilibrium E 0 , we make use of the method discussed in Castillo-Chavez et al. [7]. The general form of the theorem is stated in Appendix.

Theorem 2.3

The infection-free equilibrium point E 0 of Systems (2.18)–(2.20) is globally asymptotically stable whenever R 0 < 1 .

Proof

We will prove the global stability of E 0 = ( ω μ , 0 , 0 ) of Systems (2.18)–(2.20) by showing that Systems (2.18)–(2.20) can be written in the general form stated in Theorem A.2 (Appendix) and both the conditions A 1 and A 2 of Theorem A.2 are satisfied.

Comparing the general system to Systems (2.18)–(2.20), the functions F and G are given by

F ( X , Y ) = Λ − β N h S I − μ S ,

G ( X , Y ) = ( β N h S I − ( m u + π + γ 1 ) E , π E − ( μ + γ 2 + d N h ) I ) ,

where X = S and Y = ( E , I ) .

The infection-free equilibrium point is U 0 = ( X 0 , 0 ¯ ) , where

X 0 = Λ μ and 0 ¯ = ( 0 , 0 ) .

From the stability analysis of E 0 , we know that U 0 is locally asymptotically stable iff R 0 < 1 . Clearly, we see that G ( X , 0 ¯ ) = ( 0 , 0 ¯ ) . Now, we show that X 0 = ( Λ μ ) is globally asymptotically stable for the subsystem

(2.29) d S d t = F ( S , 0 ¯ ) = ω − μ S .

The integrating factor is e μ t , and therefore, after performing integration on equation (2.29), we obtain

S ( t ) e μ t = Λ e μ t μ + c .

As t → ∞ , we obtain

S ( t ) = ω μ ,

which is independent of c. This independency implies that X 0 = Λ μ 1 is globally asymptotically stable for the subsystem d S d t = Λ − μ S . So, the assumption A 1 is satisfied.

Now, we will show that assumption A 2 holds. First, we will find the matrix A . As per the theorem, A = D Y G ( X , Y ) at X = X 0 and Y = 0 ¯ . Now,

D Y G ( X , Y ) = − ( μ + π + γ 1 ) β N h S π − ( μ + γ 2 + d N h ) .

At X = X 0 and Y = 0 ¯ , we obtain

A = − ( μ + π + γ 1 ) β N h Λ μ π − ( μ + γ 2 + d N h ) .

Clearly, matrix A has non-negative off-diagonal elements. Hence, A is a M-matrix. Using G ^ ( X , Y ) = A Y − G ( X , Y ) , we obtain

G ^ ( X , Y ) = G 1 ^ ( X , Y ) G 2 ^ ( X , Y ) = β N h I Λ μ − S 0

Hence, G 1 ^ ( X , Y ) = β N h I ( Λ μ − S ) ≥ 0 because S ( t ) ≤ Λ μ and G 2 ^ ( X , Y ) = 0 .

Thus, both Assumptions A 1 and A 2 are satisfied, and therefore, infection-free equilibrium point E 0 is globally asymptotically stable provided R 0 < 1 .□

2.3.2 Stability analysis of E 1

The Jacobian matrix of Systems (2.18)–(2.20) at E 1 is given by

J = − ( β N h I * + μ ) 0 − β N h S β N h I * − ( μ + π + γ 1 ) β N h S * 0 π − ( μ + γ 2 + d N h ) .

The characteristic equation of the Jacobian J evaluated at E 1 is given by

(2.30) λ 3 + A 1 λ 2 + B 1 λ + C 1 = 0 ,

where

A 1 = 3 μ + π + β N h I * + γ 1 + γ 2 + d N h ,

B 1 = ( μ + β N h I * ) ( 2 μ + π + γ 1 + γ 2 + d N h ) + ( μ + π + γ 1 ) ( μ + γ 2 + d N h ) − β N h π S * ,

C 1 = ( μ + β N h I * ) ( ( μ + π + γ 1 ) ( μ + γ 2 + d N h ) − β N h π S * ) + β 2 N h 2 S * I * π .

Clearly, A 1 > 0 . By the Routh-Hurwitz criterion, all the roots of characteristic equation (2.30) are negative iff C 1 > 0 and A 1 B 1 − C 1 > 0 .

Simplifying the expression for C 1 , we obtain

C 1 = μ ( μ + π + γ 1 ) ( μ + γ 2 + d N h ) ( R 0 − 1 ) .

Therefore, C 1 > 0 iff R 0 > 1 . Hence, we conclude that the infected equilibrium point E 1 exists and remains locally asymptotically stable provided R 0 > 1 and ( A 1 B 1 − C 1 ) > 0 . In the following theorem, we summarize the aforementioned discussion on the stability of E 1 .

Theorem 2.4

A unique infected equilibrium point E 1 exists for Systems (2.18)–(2.20) and remains locally asymptotically stable provided that the following conditions are satisfied:

R 0 > 1 .
( A 1 B 1 − C 1 ) > 0 .

2.4 Bifurcation analysis

We use the method given by Castillo-Chavez and Song in [8] to do the bifurcation analysis.

Applying Theorem A.3 (appendix) to the system ( 2.18 ) – ( 2.20 ) ̲ :

In our case, we have x = ( S , E , I ) ∈ R 3 , where x 1 = S , x 2 = E , and x 3 = I . Let us consider β (transmission rate of the infection) to be the bifurcation parameter.

We know that

R 0 = β N h Λ π μ ( μ + π + γ 1 ) ( μ + γ 2 + d N h ) .

Therefore, we have

β = R 0 μ ( μ + π + γ 1 ) ( μ + γ 2 + d N h ) β N h Λ π .

Let β = β * at R 0 = 1 . So, we have

β * = μ ( μ + π + γ 1 ) ( μ + γ 2 + d N h ) β N h Λ π .

With x = ( x 1 , x 2 , x 3 ) = ( S , I , V ) Systems (2.18)–(2.20) can be written as follows:

d x 1 d t = Λ − β N h x 1 x 3 − μ x 1 = f 1 ,

d x 2 d t = β N h x 1 x 3 − ( μ + π + γ 1 ) x 2 = f 2 ,

d x 3 d t = π x 2 − ( μ + γ 2 + d N h ) x 3 = f 3 .

The infection-free equilibrium point E 0 is given by

x * = Λ μ , 0 , 0 = ( x 1 * , x 2 * , x 3 * ) .

Clearly, f ( x * , β ) = 0 , ∀ β ∈ R , where f = ( f 1 , f 2 , f 3 ) . Let D x f ( x * , β * ) denote the Jacobian matrix of the aforementioned system at the equilibrium point x * and R 0 = 1 . Now, we see that

D x f ( x * , β * ) = − μ 0 − β * N h Λ μ 0 − ( μ + π + γ 1 ) β * N h Λ μ 0 π − ( μ + γ 2 + d N h ) .

The characteristic polynomial of the aforementioned matrix D x f ( x * , β * ) is given by

(2.31) ( − μ − λ ) ( ( μ + π + γ 1 ) + λ ) ( ( μ + γ 2 + d N h ) + λ ) − β * N h Λ π μ = 0 .

Hence, we obtain the first eigenvalue of (2.31) as

λ 1 = − μ < 0 .

The other eigenvalues λ 2 , 3 of (2.31) are the solutions of the following equation:

(2.32) λ 2 + ( 2 μ + π + γ 1 + γ 2 + d N h ) λ + ( μ + π + γ 1 ) ( μ + γ 2 + d N h ) − β * N h π Λ μ = 0 .

Substituting the expression for β * in (2.32), we obtain

(2.33) λ 2 + ( 2 μ + π + γ 1 + γ 2 + d N h ) λ = 0 .

The eigenvalues of (2.33) are λ 2 = 0 and λ 3 = − ( 2 μ + π + γ 1 + γ 2 + d N h ) .

Hence, the matrix D x f ( x * , β * ) has zero as its simple eigenvalue and all other eigenvalues with negative real parts. Thus, Condition 1 of Theorem A.3 (appendix) is satisfied.

Next, for proving Condition 2, we need to find the right and left eigenvectors of the zero eigenvalue ( λ 2 ). Let us denote the right and left eigen vectors by w and v , respectively. To find w , we use ( D x f ( x * , β * ) − λ 2 I d ) w = 0 , which implies that

− μ 0 − β * N h Λ μ 0 − ( μ + π + γ 1 ) β * N h Λ μ 0 π − ( μ + γ 2 + d N h ) w 1 w 2 w 3 = 0 0 0 ,

where w = ( w 1 , w 2 , w 3 ) T . As a result, we obtain the system of simultaneous equations as follows:

(2.34) − μ w 1 − β * N h Λ μ w 3 = 0 ,

(2.35) − ( μ + π + γ 1 ) w 2 + β * N h Λ μ w 3 = 0 ,

(2.36) π w 2 − ( μ + γ 2 + d N h ) w 3 = 0 .

By choosing w 3 = μ in the aforementioned simultaneous equations (2.34)–(2.36), we obtain

w 2 = β * N h Λ ( μ + π + γ 1 ) and w 1 = − β * N h Λ μ .

Therefore, the right eigen vector of zero eigenvalue is given by

w = − β * N h Λ μ , β * N h Λ ( μ + π + γ 1 ) , μ .

Similarly, to find the left eigenvector v , we use v ( D x f ( x * , β * ) − λ 2 I d ) = 0 , which implies that

v 1 v 2 v 3 − μ 0 − β * N h Λ μ 0 − ( μ + π + γ 1 ) β * N h Λ μ 0 π − ( μ + γ 2 + d N h ) = 0 0 0 ,

where v = ( v 1 , v 2 , v 3 ) . The simultaneous equations obtained thereby are as follows:

(2.37) − μ v 1 = 0 ,

(2.38) − ( μ + π + γ 1 ) v 2 + π v 3 = 0 ,

(2.39) − β * N h Λ μ v 1 + β * N h Λ μ v 2 − ( μ + γ 2 + d N h ) v 3 = 0 .

Therefore, solving the aforementioned simultaneous equations (2.37)–(2.39), we obtain v 1 = 0 .

By choosing v 2 = 1 , we obtain

v 3 = β * N h Λ ( μ ( μ + γ 2 + d N h ) ) .

Hence, the left eigen vector is given by

v = 0 , 1 , β * N h Λ ( μ ( μ + γ 2 + d N h ) ) .

Now, we need to find a and b . As per Theorem 2.7, a and b are given by

a = ∑ k , i , j = 1 3 v k u i u j ∂ 2 f k ∂ x i ∂ x j ( x * , β * ) ,

b = ∑ k , i = 1 3 v k u i ∂ 2 f k ∂ x i ∂ β ( x * , β * ) .

Expanding the summation in the expression for a , it reduces to

a = w 1 w 3 ∂ 2 f 2 ∂ x 1 ∂ x 3 + w 3 w 1 ∂ 2 f 2 ∂ x 3 ∂ x 1 ,

where partial derivatives are found at ( x * , β * ) . Now,

∂ 2 f 2 ∂ x 1 ∂ x 3 ( x * , β * ) = β * N h and ∂ 2 f 2 ∂ x 3 ∂ x 1 ( x * , β * ) = β * N h .

Substituting these partial derivatives along with w in the expression of a , we obtain

a = − 2 β * 2 N h 2 Λ < 0 .

Next, expanding the summation in the expression for b , we obtain

b = v 2 w 2 ∂ 2 f 2 ∂ x 1 ∂ β ( x * , β * ) + v 2 w 3 ∂ 2 f 2 ∂ x 3 ∂ β ( x * , β * ) .

Now,

∂ 2 f 2 ∂ x 3 ∂ β ( x * , β * ) = N h λ μ ,

∂ 2 f 2 ∂ x 1 ∂ β ( x * , β * ) = 0 .

Substituting in the expression of b , we obtain

b = N h Λ > 0 .

Hence, a < 0 and b > 0 . We note that Condition (iv) of Theorem A.3 (appendix) is satisfied. Hence, we conclude that the system undergoes bifurcation at β = β * , implying R 0 = 1 .

Thus, we conclude that when R 0 < 1 , there exists a unique infection-free equilibrium, which is globally asymptotically stable, and negative infected equilibrium which is unstable. Since negative values of population are not practical, we ignore it in this case. Furthermore, as R 0 crosses unity from below, the infection-free equilibrium point loses its stable nature and become unstable, the bifurcation point being at β = β * implying R 0 = 1 and there appears a positive locally asymptotically stable infected equilibrium point. There is an exchange of stability between infection-free equilibrium and infected equilibrium at R 0 = 1 . Hence, a forward bifurcation (trans-critical bifurcation) takes place at the break point β = β * . We summarize the aforementioned discussion on bifurcation by the following theorem.

Theorem 2.5

As R 0 crosses unity, the infection-free equilibrium changes its stability from stable to unstable and there exists a locally asymptotically infected equilibrium when R 0 > 1 , i.e., direction of bifurcation is forward (transcritical) at R 0 > 1 .

From the expression of R 0 , the condition R 0 = 1 can be rewritten in terms of the total viral load N h at within-host scale. From equation (2.26), we see that R 0 = 1 implies N h = μ ( μ + π + γ 1 ) ( μ + γ 2 ) β π Λ − μ d ( μ + π + γ 1 ) . Therefore, we infer from this that, there is a critical value of the viral load at within-host scale depending on which the stability of the reduced multi-scale model (2.18)–(2.20) changes.

2.5 Sensitivity and elasticity

The basic reproduction number denoted by R 0 is one of the most important quantities in any infectious disease models. The expression for R 0 for the reduced multi-scale model (2.18)–(2.20) is given by

R 0 = β N h Λ π μ ( μ + π + γ 1 ) ( μ + γ 2 + d N h ) .

To determine the best control measures, knowledge of the relative importance of the different factors responsible for transmission is useful. Initially, disease transmission is related to R 0 and sensitivity predicts which parameters have a high impact on R 0 . The sensitivity index of R 0 with respect to a parameter μ is ∂ R 0 ∂ μ : another measure is the elasticity index (normalized sensitivity index) that measures the relative change of R 0 with respect to μ , denoted by ϕ μ R 0 , and defined as

ϕ μ R 0 = ∂ R 0 ∂ μ μ R 0 .

The sign of the elasticity index tells whether R 0 increases (positive sign) or decreases (negative sign) with the parameter; however, the magnitude determines the relative importance of the parameter [5,42]. These indices can guide control by indicating the most important parameters to target, although feasibility and cost play a role in practical control strategy. If R 0 is known explicitly, then the elasticity index for each parameter can be computed explicitly and evaluated for a given set of parameters.

The elasticity index of R 0 with the model parameters is given by

ϕ β R 0 = ∂ R 0 ∂ β β R 0 = 1 ,

ϕ Λ R 0 = ∂ R 0 ∂ Λ Λ R 0 = 1 ,

ϕ π R 0 = ∂ R 0 ∂ π π R 0 = 1 ,

ϕ μ R 0 = ∂ R 0 ∂ μ μ R 0 = − β N h Λ π ( 3 μ 2 + 2 μ ( π + γ 1 + γ 2 + d N h ) + ( γ 2 + d ) ( π + γ 1 ) ) ( μ 3 + μ 2 ( π + γ 1 + γ 2 + d N h ) + μ ( γ 2 + d ) ( π + γ 1 ) ) μ R 0 ,

ϕ γ 1 R 0 = ∂ R 0 ∂ γ 1 γ 1 R 0 = − γ 1 ( μ + π + γ 1 ) ,

ϕ γ 2 R 0 = ∂ R 0 ∂ γ 2 γ 2 R 0 = − γ 2 ( μ + γ 1 + d N h ) ,

ϕ N h R 0 = ∂ R 0 ∂ N h N h R 0 = − μ + γ 2 ( μ + γ 2 + d N h ) ,

ϕ d R 0 = ∂ R 0 ∂ d d R 0 = − d N h ( μ + γ 2 + d N h ) .

The elasticity indices calculated using the parameter values from Table 3 are given in Table 2. The elastic index of parameters β , π , Λ , and N h are positive, and the remaining are negative. This implies that the increase in the values of these parameters increases R 0 , whereas increase in the values of parameters μ , γ 1 , and γ 2 decreases R 0 . For parameter β , ϕ β R 0 = 1 implies an increase (decrease) of β by y % increases (decreases) R 0 by the same percentage. From Table 2, we see that the basic reproduction number is most sensitive to the parameters β , π , and Λ . The implication of this is that an increase in the transmission rate increases the spread of the disease in the community.

Table 2

Elasticity indices of R 0

Parameters	Elastic index	Value
β	ϕ β R 0	1
π	ϕ π R 0	1
Λ	ϕ Λ R 0	1
μ	ϕ μ R 0	‒ 0.2785
γ 1	ϕ γ 1 R 0	‒ 0.3196
γ 2	ϕ γ 2 R 0	‒ 0.00082
N h	ϕ N h R 0	0.0018
d	ϕ d R 0	− 0.99

2.6 Numerical illustrations

2.6.1 Numerical illustrations of the stability of equilibrium points

Now, we numerically illustrate the stability of the equilibrium points admitted by Systems (2.18)–(2.20). The simulation is done using MatLab software and ode solver ode45 is used to solve the system of equations. The parameter values of the within-host sub-model used in simulation are taken from [10]. These within-host parameter values are used in calculating the area under the viral load curve N h . The parameter values of the reduced multi-scale model (2.18)–(2.20) are taken from [34,39]. All the parameter values are listed in Table 3. The initial values of the variables are given in Table 4. These values are used in illustrating the stability of the equilibrium points admitted by Systems (2.18)–(2.20). We also illustrate the influence of key within-host scale sub-model parameters on the between-host scale sub-model variables.

Table 3

Parameter values

Symbols	Values	Source
ω	2	[10]
k	0.05	[10]
μ c	0.1	[10]
μ v	0.1	[10]
α	0.24	[28]
x	0.795	[10]
y	0.56	[10]
Λ	μ N ( 0 ) = 71.3	[39]
β	0.0115	[39]
μ	0.062	[34]
π	0.09	[39]
d	0.0018	[39]
γ 1 , γ 2	0.05, 0.0714	[39]

Table 4

Initial values of the variables

Variable	Initial values	Source
U ( s )	3.2 × 1 0 5	[10]
U * ( s )	0	[10]
V ( s )	5.2	[10]
S ( t )	1,000	Assumed
E ( t )	100	Assumed
I ( t )	50	Assumed

Taking d 1 = 0 and d 2 = 30 and the other parameter values from Table 3, the area under the viral load curve N h is found to be 3.3759 × 1 0 4 . Here, in our study, we have taken d 2 − d 1 = 30 days as the duration of host infectivity. The value of basic reproduction number R 0 with β = 0.00115 and μ = 0.42 , and other parameters from Table 3 are calculated to be R 0 = 0.84 . Since R 0 < 1 , by Theorem 2.3, the infection-free equilibrium point, E 0 = ( 169.76 , 0 , 0 ) is globally asymptotically stable for Systems (2.18)–(2.20). The plots in Figures 1 and 2 for different initial conditions depict the local and global asymptotic stability of the infection-free equilibrium E 0 of Systems (2.18)–(2.20).

$Figure 1 Local asymptotic stability of E 0 {E}_{0} of Systems (2.18)–(2.20).$

Figure 1

Local asymptotic stability of E 0 of Systems (2.18)–(2.20).

$Figure 2 Global asymptotic stability of E 0 {E}_{0} of Systems (2.18)–(2.20).$

Figure 2

Global asymptotic stability of E 0 of Systems (2.18)–(2.20).

We know that the infected equilibrium E 1 exists only if R 0 > 1 and it is also locally asymptotically stable whenever R 0 > 1 . For the parameter values in Table 3, the value of R 0 was calculated to be 135.7936 and E 1 to be ( 8.4687 , 316.8 , and 165.03 ) . Figure 3 demonstrates that E 1 is locally asymptotically stable whenever R 0 > 1 . From Section 2.4, we know that the equilibrium point exchanges their stabilities depending on the value of R 0 . In Figure 4, we show the occurrence of the trans-critical bifurcation at R 0 = 1 .

$Figure 3 Local asymptotic stability of E 1 {E}_{1} of Systems (2.18)–(2.20) starting from the initial state ( S 0 , E 0 , I 0 ) = ( 1,000 , 100 , 50 ) . \left({S}_{0},{E}_{0},{I}_{0})=\left(\hspace{0.1em}\text{1,000}\hspace{0.1em},\hspace{0.33em}100,50).$

Figure 3

Local asymptotic stability of E 1 of Systems (2.18)–(2.20) starting from the initial state ( S 0 , E 0 , I 0 ) = ( 1,000 , 100 , 50 ) .

$Figure 4 Trans-critical bifurcation exhibited by Systems (2.18)–(2.20) at R 0 = 1 . {R}_{0}=1. The change in the stability of the equilibria with variation in R 0 {R}_{0} can be observed.$

Figure 4

Trans-critical bifurcation exhibited by Systems (2.18)–(2.20) at R 0 = 1 . The change in the stability of the equilibria with variation in R 0 can be observed.

2.6.2 Heat plots

Here, we vary two parameters of the model (2.18)–(2.20) at a time in a certain interval and plot the value of R 0 as heat plots. Heat plot has two different colors: one corresponding to the region with R 0 < 1 and the other corresponding to R 0 > 1 . This plot helps to find the combination of parameter values for which R 0 < 1 and R 0 > 1 . The blue color region in these plots corresponds to the region where R 0 < 1 , and therefore, from Theorem 2.3, the infection-free equilibrium is globally stable in this region. The other region with green color corresponds to the region where R 0 > 1 , and in this region, the infection-free equilibrium is unstable and there exists a unique infected equilibrium point whose stability is determined using Theorem 2.4. The stability of infected equilibria is determined using Theorem 2.4.

Parameters μ and d ̲

Heat plot varying μ in the interval ( 0.5 , 0.94 ) and d in the interval ( 0.001 , 0.44 ) is given in Figure 5. When the values of the parameters μ and d fall in the blue colored region, the infection-free equilibrium point E 0 of Systems (2.18)–(2.20) is globally asymptotically stable, and in the green colored region, it is unstable.

$Figure 5 Heat plots for parameters μ \mu and d d .$

Figure 5

Heat plots for parameters μ and d .

Parameters π and d ̲

Heat plot varying π in the interval ( 0.0005 , 0.005 ) and d in the interval ( 0.05 , 0.5 ) is given in Figure 6. The infection-free equilibrium point E 0 of Systems (2.18)–(2.20) is globally asymptotically stable, in the blue colored region, and unstable in the green colored region.

$Figure 6 Heat plots for parameters π \pi and d d .$

Figure 6

Heat plots for parameters π and d .

Parameters μ and π ̲

Heat plot varying π in the interval ( 0.004 , 0.049 ) and μ in the interval ( 1.05 , 1.5 ) is given in Figure 7 and similarly varying β in the interval ( 0.0005 , 0.005 ) and d in the interval ( 0.05 , 0.5 ) is given in Figure 8.

$Figure 7 Heat plots for parameters μ \mu and π \pi .$

Figure 7

Heat plots for parameters μ and π .

$Figure 8 Heat plots for parameters β \beta and d d .$

Figure 8

Heat plots for parameters β and d .

2.6.3 Influence of within-host sub-model parameters on the between-host sub-model

In this section, we study the influence of the within-host sub-model parameters on the between-host sub-model variables. The reduced multi-scale model (2.18)–(2.20) is categorized as an NMSM according to the categorization of multi-scale model infectious disease systems [15]. Therefore, the multi-scale model (2.18)–(2.20) is unidirectionally coupled. In this only the within-host scale sub-model influences the between-host scale sub-model without any reciprocal feedback. Here, we illustrate the influence of the key within-host scale sub-model parameters such as α , y , and x on the between-host scale sub-model variable I . At the within-host scale sub-model, α , y , and x describe the production of virus by infected cells, the clearance of the infected cells by the immune system, and the clearance rate of free virus particles by the immune system. The parameter values for the simulation are taken from Table 3.

In Figure 9, the effect of variation of burst rate of virus on infected population is plotted. The infected population I ( t ) is plotted for three different values of burst rate α = 0.24 , α = 0.5 , and α = 0.7 . The values of N h will change based on the values of α . This change in the value of N h will also affect the between-host dynamics. The values of N h for α = 0.24 , α = 0.5 , and α = 0.7 were calculated to be 3.3759 × 1 0 4 , 3.6921 × 1 0 4 , and 3.8127 × 1 0 4 , respectively. Figure 9 shows that the between-host scale variable I ( t ) is influenced by the within-host scale parameter α . We see from Figure 9 that as the infected cell burst size increases, transmission of the disease in the community also increases. Therefore, antiviral drugs such as remdesivir, arbidol, and hydroxychloroquine (HCQ) that reduce the average infected cell viral production rate at within-host scale will possibly reduce the transmission of COVID-19 disease at between-host scale.

Figure 9

Effect of variation of burst rate of virus on infected population.

In Figure 10, the effect of variation of rate of clearance of infected cells by the immune system on infected population is plotted. The infected population I ( t ) is plotted for three different values x = 0.5 , x = 0.795 , and x = 0.85 . The values of N h were calculated to be 3.5149 × 1 0 4 , 3.3759 × 1 0 4 , and 2.8767 × 1 0 4 , respectively. These numerical illustrations show that the between-host scale variable I ( t ) is influenced by the within-host scale parameter x . We see from Figure 10 that as rate of clearance of infected cells by the immune system increases, the infection at population level decreases. Therefore, drugs that kill infected cells could have community-level benefits of reducing COVID-19 transmission.

Figure 10

Effect of variation of rate of clearance of infected cells by the immune system on infected population.

In Figure 11, the effect of variation of rate of clearance of virus particles by the immune system on infected population is plotted. The infected population I ( t ) is plotted for three different values y = 0.56 , y = 0.7 , and y = 0.8 . The values of N h for y = 0.56 , y = 0.7 , and y = 0.8 were calculated to be 3.3612 × 1 0 4 , 3.211 × 1 0 4 , and 3.1227 × 1 0 4 , respectively. These numerical illustrations show that the between-host scale variable I ( t ) is influenced by the within-host scale parameter y . We see from Figure 11 that as rate of clearance of virus by the immune system increases, the infection in the community decreases. Therefore, in addition to the benefits at individual level, treatments that increase the clearance rate of free virus particles in an infected individual could also have potential benefits in reducing the transmission at between-host level.

Figure 11

Effect of variation of rate of clearance of virus by the immune system on infected population.

3 Comparative effectiveness of health interventions

In this section, in similar lines to the study done in [16,36], we extend the work on the multi-scale model discussed in Section 2 by including health care interventions. In the context of infectious diseases, disease dynamics can cause a large difference between the performance of a health intervention at the individual level (within-host scale), which is easy to determine, and its performance at the population level (between-host scale), which is difficult to determine [16]. As a result, in situations where the effectiveness of a health intervention cannot be determined, health interventions with proven effectiveness may be recommended over those with potentially higher comparable effectiveness [16]. We use three health care interventions that are as follows:

1. Antiviral drugs:

Drugs such as remdesivir inhibit RNA-dependent RNA polymerase, and drugs lopinavir and ritonavir inhibit the viral protease, thereby reducing viral replication [41]. Considering that these antiviral drugs are administered, the burst rate of the within-host scale sub-model α now gets modified to α ( 1 − ε ) , where ε is the efficacy of antiviral drugs and 0 < ε < 1 .
HCQ acts by preventing the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) virus from binding to cell membranes and also chloroquine/HCQ could inhibit viral entry by acting as inhibitors of the biosynthesis of sialic acids, critical actors of virus–cell ligand recognition [38]. The administration of this drug decreases the infection rate k of susceptible cells. This parameter gets modified to become k ( 1 − γ ) , where γ is the efficacy of these drugs and 0 < γ < 1 .

2. Immunomodulators:

Interferons are broad-spectrum antivirals, exhibiting both direct inhibitory effect on viral replication and supporting an immune response to clear virus infection [43]. Due to the administration of immunomodulators such as INF, immunity power of an individual increases as a result of which the clearance rate of infected cells and virus particles by immune response increases. Therefore, the parameters x and y now get modified to x ( 1 + δ ) and y ( 1 + δ ) , where δ is the efficacy of immunomodulators and 0 < δ < 1 .

3. Generalized social distancing:

This intervention involves introducing measures such as restriction of mass gatherings, wearing of mask, and temporarily closing schools. Assuming that generalized social distancing is introduced to control COVID-19 epidemic, then the rate of contact with community β obtains modified to β ( 1 − ρ ) , where ρ is the efficacy of generalized social distancing and 0 < ρ < 1 .

Considering all the modifications in the parameters, the multi-scale model incorporating the effects of all the aforementioned health interventions becomes

(3.1) d S d t = Λ − β ( 1 − ρ ) N m S ( t ) I ( t ) − μ S ( t ) ,

(3.2) d E d t = β N m S ( t ) I ( t ) − ( μ + π + γ 1 ) E ( t ) ,

(3.3) d I d t = π E ( t ) − ( μ + γ 2 ) I ( t ) − d N m I ( t ) ,

where N m is modified N h given by

(3.4) N m = α ( 1 − ε ) ∫ d 1 d 2 U * d s y ( 1 + δ ) + μ v .

The effective reproduction number after incorporating the health interventions is given by

R E = β ( 1 − ρ ) N m Λ π μ ( μ + π + γ 1 ) ( μ + γ 2 + d N m ) .

Here, the comparative effectiveness of the aforementioned three health interventions are evaluated using R E as indicators of intervention effectiveness. In order to evaluate these, we calculate the percentage reduction of ℛ 0 for single and multiple combination of the interventions at different efficacy levels such as (a) low efficacy of 0.3, (b) medium efficacy of 0.6, and (c) high efficacy of 0.9. The different efficacy levels of the health interventions are chosen from [36].

Percentage reduction of ℛ 0 is given by

ℛ 0 − ℛ E j ℛ 0 × 100 ,

where R E j means the effective reproductive number when intervention/combination of interventions with efficacy/efficacies j .

We now consider eight different combinations of these three health interventions corresponding to efficacy values 0.3, 0.6, and 0.9 obtained using the effective reproductive number R E as the indicator of intervention effectiveness. Then, for each efficacy level, we rank the percentage reductions on ℛ 0 in ascending order from 1 to 8 corresponding to the different combinations of three health interventions considered in this study. The comparative effectiveness is calculated and measured on a scale from 1 to 8, with 1 denoting the lowest comparative effectiveness and 8 denoting the highest comparative effectiveness. In Table 5, the abbreviations (a) CEL stands for “comparative effectiveness at low efficacy,” which is 0.3; (b) CEM stands for “comparative effectiveness at medium efficacy,” which is 0.6; and (c) CEH stands for “comparative effectiveness at high efficacy” which is 0.9.

Table 5

Comparative effectiveness for ℛ 0

No.	Indicator	% Age	CEL	% Age	CEM	% Age	CEH
1	ℛ 0	0	1	0	1	0	1
2	ℛ E ρ	30	5	60	5	90	5
3	ℛ E δ	0.106	3	0.24	2	0.45	4
4	ℛ E ε	0.075	2	0.26	3	0.38	3
5	ℛ E ρ δ	30.07	7	60.12	7	90.45	7
6	ℛ E ρ ε	30.05	6	60.1	6	90.23	6
7	ℛ E ε δ	0.22	4	0.35	4	0.67	2
8	ℛ E ρ δ ε	30.16	8	60.34	8	90.5	8

From Table 5, we deduce the following results regarding comparative effectiveness of three different interventions considered:

When a single intervention strategy is implemented, intervention involving introducing measures such as restriction of mass gatherings, wearing of mask, and temporarily closing schools show a significant decrease in R 0 compared to the implementation of antiviral drugs and immunomodulators at all efficacy levels.
When considering two interventions, we observe that the generalized social distancing along with immunomodulators that boost the immune response would be highly effective in limiting the spread of infection in the community.
A combined strategy involving treatment with antiviral drugs, immunomodulators, and generalized social distancing seems to perform the best among all the combinations considered at all efficacy levels.

From this section, we conclude that a combined strategy is the best strategy to contain the spread of infection. However, the implementation of a combined strategy may not always be cost-effective and feasible. In a single intervention strategy, generalized social distancing is shown to lower the value of R 0 better than the individual use of antiviral drugs and immunomodulators. Therefore, in a situation where resources are limited and costly, interventions such as restricting mass gatherings, wearing masks, and temporarily closing schools would be highly effective in containing infection and also incur less cost.

4 Discussions and conclusions

COVID-19 is a contagious respiratory and vascular disease caused by SARS-CoV-2. The emergence of the disease has posed a major challenge to health authorities and governments around the world. The COVID -19 has spread globally and is one of the largest pandemics the world has ever seen. Mathematical models have proved to provide useful information about the dynamics of the infectious diseases in the past. To understand the dynamics of the COVID-19 disease, several models have been developed and studied [6,9,12,27,35,39,40,44,47,48]. Although these models have provided useful information about the disease transmission, disease modelers see the need to integrate or link the models at different scales (within-host and between-host) to obtain a better understanding of the dynamics of infectious disease spread.

In this study, we develop an NMSM for COVID-19 disease that integrates the within-host scale and the between-host scale sub-models. The transmission rate and COVID-19-induced death rate at between-host scale are assumed to be a linear function of viral load. Because of the difficulties of working at two different time scales, we approximate individual-level host infectiousness by some surrogate measurable quantity called area under viral load and reduce the multi-scale model at two different times scales to a model with area under the viral load curve acting as a proxy of individual level host infectiousness.

Initially, the well-posedness of the reduced multi-scale model is discussed followed by the stability analysis of the equilibrium points admitted by the reduced multi-scale model. The infection-free equilibrium point of the model is found to remain globally asymptotically stable whenever the value of basic reproduction number is less than unity. As the value of basic reproduction number crossed unity, a unique infected equilibrium point exists and remains asymptotically stable if A 1 B 1 > C 1 (Theorem 2.4). The model is shown to undergo a forward (trans-critical) bifurcation at R 0 = 1 . To predict the sensitivity of the model parameters on R 0 , elastic index that measures the relative change of R 0 with respect to parameters is calculated for each parameters in the definition of R 0 . The parameters β , π , and Λ were found to be most sensitive toward R 0 . The theoretical results are supported with numerical illustrations. To separate out the region of stability and instability of the equilibrium points, two-parameter heat plot is done by varying the parameter values in certain range. The influence of the key within-host scale sub-model parameters such as α , y , and x on between-host scale are also numerically illustrated. It is found that the spread of infection in a community is influenced by the production of virus particles by an infected cells ( α ), clearance rate of infected cell ( x ) , and clearance rate of virus particles by immune system ( y ) .

We also use the reduced multi-scale model developed to study the comparative effectiveness of the three health interventions (antiviral drugs, immunomodulators, and generalized social distancing) for COVID-19 viral infection using R E as the indicator of intervention effectiveness. The result suggested that a combined strategy involving treatment with antiviral drugs, immunomodulators, and generalized social distancing would be the best strategy to limit the spread of infection in the community. However, implementing a combined strategy may not be always cost-effective or feasible. In a single intervention strategy, general social distancing has been shown to lower the value of R 0 better than the individual use of antiviral drugs and immunomodulators. Therefore, in a situation where resources are limited and costly, interventions such as restricting mass gatherings, wearing masks, and temporarily closing schools would be highly effective in containing infection and also more cost-effective.

The study performed at multi-scale level in this study highlights the dependence of between-host processes on within-host processes. With the increasing recognition of the importance of multi-scale modeling in disease dynamics, we believe that our study contributes to the growing body of knowledge on multi-scale modeling of COVID-19 disease. Two categories of people benefit from the multi-scale study: disease modelers and public health officials. Disease modelers benefit from the fact that this study provides a different approach to modeling acute viral infections. For public health planners, this study provides a model-based approach to evaluating the effectiveness of public health interventions. We also believe that the results presented in this study will help physicians, doctors, and researchers in making informed decisions about COVID-19 disease prevention and treatment interventions.

Funding information: This research received no specific grant from any funding agency, commercial, or nonprofit sectors.
Author contributions: BC conceived the idea and mathematical model, conducted the analysis, numerical simulations, and conceptualized the entire article. KKVD contributed to the conceptualization of the mathematical model.
Conflict of interest: The authors have no conflicts of interest to disclose.
Data availability statement: Not applicable.
Ethical approval: This research did not required ethical approval.
Dedication: The authors respectfully dedicate this paper to the Founder Chancellor of SSSIHL, Bhagawan Sri Sathya Sai Baba. Additionally, the first author dedicates this work to his father, late Purna Chhetri, and the senior author dedicates it to his loving elder brother, D.A.C. Prakash, who still lives in his heart.
Code availability: Not applicable.

Appendix

Theorem A.1

Let D denote the domain:

∣ t − t 0 ∣ ≤ a , ∣ ∣ x − x 0 ∣ ∣ ≤ b , x = ( x 1 , x 2 , … , x n ) , x 0 = ( x 10 , … , x n 0 ) ,

and suppose that f ( t , x ) satisfies the Lipschitz condition:

(A1) ∣ ∣ f ( t , x 2 ) − f ( t , x 1 ) ∣ ∣ ≤ k ∣ ∣ x 2 − x 1 ∣ ∣ ,

and whenever the pairs ( t , x 1 ) and ( t , x 2 ) belong to the domain D, where, k is used to represent a positive constant. Then, there exist a constant δ > 0 such that there exists a unique (exactly one) continuous vector solution x ( t ) of System (2.21) in the interval ∣ t − t 0 ∣ ≤ δ . It is important to note that Condition (2.22) is satisfied by requirement that

∂ f i ∂ x j , i , j = 1 , 2 , … , n ,

be continuous and bounded in the domain D.

Theorem A.2

Consider the following general system:

(A2) d X d t = F ( X , Y ) , d Y d t = G ( X , Y ) ,

where X denotes the uninfected population compartments and Y denotes the infected population compartments, including latent and infectious. Here, the function G is such that it satisfies G ( X , 0 ) = 0 . Let U 0 = ( X 0 , 0 ¯ ) denote the equilibrium point of the aforementioned general system.

If the following two conditions are satisfied, then the infection-free equilibrium point U 0 is globally asymptotically stable for the aforementioned general system provided R 0 < 1

A 1 : For the subsystem d X d t = F ( X , 0 ) , X 0 is globally asymptotically stable.

A 2 : The function G = G ( X , Y ) can be written as G ( X , Y ) = A Y − G ^ ( X , Y ) , where G j ^ ( X , Y ) ≥ 0 ∀ ( X , Y ) in the biologically feasible region Ω for j = 1 , 2 and A = D Y G ( X , Y ) at ( X 0 , 0 ¯ ) is a M-matrix (matrix with non-negative off diagonal element).

Theorem A.3

Consider a system

d X d t = f ( X , ϕ ) ,

where X ∈ R n , ϕ ∈ R is the bifurcation parameter and f : R n × R → R n , where f ∈ C 2 ( R n , R ) . Let 0 ¯ be the equilibrium point of the system such that f ( 0 ¯ , ϕ ) = 0 ¯ , ∀ ϕ ∈ R . Let the following conditions hold:

For the matrix A = D X f ( 0 ¯ , 0 ) , zero is the simple eigenvalue, and all other eigenvalues have negative real parts.
Corresponding to zero eigenvalue, matrix A has non-negative right eigenvector, denoted as u and non-negative left eigenvectors, denoted as v.

Let f k be the k th component of f. Let a and b be defined as follows:

a = ∑ k , i , j = 1 n v k w i w j ∂ 2 f k ∂ x i ∂ x j ( 0 ¯ , 0 ) ,

b = ∑ k , i = 1 n v k w i ∂ 2 f k ∂ x i ∂ ϕ ( 0 ¯ , 0 ) .

Then, local dynamics of the system near the equilibrium point 0 ¯ is totally determined by the signs of a and b. Here are the following conclusions:

If a > 0 and b > 0 , then whenever ϕ < 0 with ∣ ϕ ∣ ≪ 1 , the equilibrium 0 ¯ is locally asymptotically stable, and moreover, there exists a positive unstable equilibrium. However, when 0 < ϕ ≪ 1 , 0 ¯ is an unstable equilibrium, and there exists a negative and locally asymptotically stable equilibrium.
If a < 0 , b < 0 , then whenever ϕ < 0 with ∣ ϕ ∣ ≪ 1 , 0 ¯ is an unstable equilibrium, whereas if 0 < ϕ ≪ 1 , 0 ¯ is locally asymptotically stable equilibrium and there exists a positive unstable equilibrium.
If a > 0 , b < 0 , then whenever ϕ < 0 with ∣ ϕ ∣ ≪ 1 , 0 ¯ is an unstable equilibrium, and there exists a locally asymptotically stable negative equilibrium. However, if 0 < ϕ ≪ 1 , 0 ¯ is stable, and there appears a positive unstable equilibrium.
If a < 0 , b > 0 , then whenever ϕ changes its value from negative to positive, the equilibrium 0 ¯ changes its stability from stable to unstable. Correspondingly, a negative equilibrium, unstable in nature, becomes positive and locally asymptotically stable.

References

[1] Aguiar, M., Anam, V., Cusimano, N., Knopoff, D., & Stollenwerk, N. (2022). Understanding COVID-19 epidemics: a multi-scale modeling approach. In: Predicting Pandemics in a Globally Connected World, Volume 1: Toward a Multiscale, Multidisciplinary Framework through Modeling and Simulation (pp. 11–42). Springer. 10.1007/978-3-030-96562-4_2Search in Google Scholar

[2] Almocera, A. E. S., Hernandez-Vargas, E. A., & Nguyen, V. K. (2018). Multiscale model within-host and between-host for viral infectious diseases. Journal of Mathematical Biology, 77(4), 1035–1057. 10.1007/s00285-018-1241-ySearch in Google Scholar PubMed

[3] Alqahtani, R. T. (2021). Mathematical model of sir epidemic system (COVID-19) with fractional derivative: stability and numerical analysis. Advances in Difference Equations, 2021(1), 1–16. 10.1186/s13662-020-03192-wSearch in Google Scholar PubMed PubMed Central

[4] Bellomo, N., Burini, D., & Outada, N. (2022). Multiscale models of COVID-19 with mutations and variants. Networks and Heterogeneous Media, 17(3), 293–310. 10.3934/nhm.2022008Search in Google Scholar

[5] Berhe, H. W., Makinde, O. D., & Theuri, D. M. (2019). Parameter estimation and sensitivity analysis of dysentery diarrhea epidemic model. Journal of Applied Mathematics, 2019, 1–13. 10.1155/2019/8465747Search in Google Scholar

[6] Biswas, S. K., Ghosh, J. K., Sarkar, S., & Ghosh, U. (2020). Covid-19 pandemic in India: a mathematical model study. Nonlinear Dynamics, 102(1), 537–553. 10.1007/s11071-020-05958-zSearch in Google Scholar PubMed PubMed Central

[7] Carlos Castillo-Chavez, Z. F., & Huang, W. (2002). On the computation of reproduction number and its role in global stability. Institute for Mathematics and Its Applications, 125(2), 229–250. 10.1007/978-1-4757-3667-0_13Search in Google Scholar

[8] Castillo-Chavez, C., & Song, B. (2004). Dynamical models of tuberculosis and their applications. Mathematical Biosciences and Engineering, 1(2), 361–404. 10.3934/mbe.2004.1.361Search in Google Scholar PubMed

[9] Chen, T.-M., Rui, J., Wang, Q.-P., Zhao, Z.-Y., Cui, J.-A., & Yin, L. (2020). A mathematical model for simulating the phase-based transmissibility of a novel coronavirus. Infectious Diseases of Poverty, 9(1), 1–8. 10.1186/s40249-020-00640-3Search in Google Scholar PubMed PubMed Central

[10] Chhetri, B., Bhagat, V. M., Vamsi, D., Ananth, V., Mandale, R., Muthusamy, S., Sanjeevi, C. B., et al. (2021). Within-host mathematical modeling on crucial inflammatory mediators and drug interventions in covid-19 identifies combination therapy to be most effective and optimal. Alexandria Engineering Journal, 60(2), 2491–2512. 10.1016/j.aej.2020.12.011Search in Google Scholar

[11] Coombs, D., Gilchrist, M. A., & Ball, C. L. (2007). Evaluating the importance of within-and between-host selection pressures on the evolution of chronic pathogens. Theoretical Population Biology, 72(4), 576–591. 10.1016/j.tpb.2007.08.005Search in Google Scholar PubMed

[12] Dashtbali, M., & Mirzaie, M. (2021). A compartmental model that predicts the effect of social distancing and vaccination on controlling COVID-19. Scientific Reports, 11(1), 1–1110.1038/s41598-021-86873-0Search in Google Scholar PubMed PubMed Central

[13] Diekmann, O., Heesterbeek, J., & Roberts, M. G. (2010). The construction of next-generation matrices for compartmental epidemic models. Journal of the Royal Society Interface, 7(47), 873–885. 10.1098/rsif.2009.0386Search in Google Scholar PubMed PubMed Central

[14] Feng, Z., Velasco-Hernandez, J., Tapia-Santos, B., & Leite, M. C. A. (2012). A model for coupling within-host and between-host dynamics in an infectious disease. Nonlinear Dynamics, 68(3), 401–411. 10.1007/s11071-011-0291-0Search in Google Scholar

[15] Garira, W. (2017). A complete categorization of multiscale models of infectious disease systems. Journal of Biological Dynamics, 11(1), 378–435. 10.1080/17513758.2017.1367849Search in Google Scholar PubMed

[16] Garira, W., & Mathebula, D. (2020). Development and application of multiscale models of acute viral infections in intervention research. Mathematical Methods in the Applied Sciences, 43(6), 3280–3306. 10.1002/mma.6119Search in Google Scholar

[17] Gilchrist, M. A., & Sasaki, A. (2002). Modeling host-parasite coevolution: a nested approach based on mechanistic models. Journal of Theoretical Biology, 218(3), 289–308. 10.1006/jtbi.2002.3076Search in Google Scholar PubMed

[18] Guo, D., Li, K. C., Peters, T. R., Snively, B. M., Poehling, K. A., & Zhou, X. (2015). Multi-scale modeling for the transmission of influenza and the evaluation of interventions toward it. Scientific Reports, 5(1), 1–9. 10.1038/srep08980Search in Google Scholar PubMed PubMed Central

[19] Gutierrez, J. B., Galinski, M. R., Cantrell, S., & Voit, E. O. (2015). From within host dynamics to the epidemiology of infectious disease: scientific overview and challenges. Mathematical Biosciences, 270, 143–155. 10.1016/j.mbs.2015.10.002Search in Google Scholar PubMed PubMed Central

[20] Hadjichrysanthou, C., Cauët, E., Lawrence, E., Vegvari, C., De Wolf, F., and Anderson, R. M. (2016). Understanding the within-host dynamics of influenza a virus: from theory to clinical implications. Journal of The Royal Society Interface, 13(119), 20160289. 10.1098/rsif.2016.0289Search in Google Scholar PubMed PubMed Central

[21] Handel, A., & Rohani, P. (2015). Crossing the scale from within-host infection dynamics to between-host transmission fitness: a discussion of current assumptions and knowledge. Philosophical Transactions of the Royal Society B: Biological Sciences, 370(1675), 20140302. 10.1098/rstb.2014.0302Search in Google Scholar PubMed PubMed Central

[22] Heldt, F. S., Frensing, T., Pflugmacher, A., Gröpler, R., Peschel, B., & Reichl, U. (2013). Multiscale modeling of influenza a virus infection supports the development of direct-acting antivirals. PLoS Computational Biology, 9(11), e1003372. 10.1371/journal.pcbi.1003372Search in Google Scholar PubMed PubMed Central

[23] Jeffery, G. M., & Eyles, D. E. (1955). Infectivity to mosquitoes of plasmodium falciparum as related to gametocyte density and duration of infection1. The American Journal of Tropical Medicine and Hygiene, 4(5), 781–789. 10.4269/ajtmh.1955.4.781Search in Google Scholar PubMed

[24] Kaplan, J. E., Khabbaz, R. F., Murphy, E. L., Hermansen, S., Roberts, C., Lal, R., …, Schreiber, G. B. (1996). Male-to-female transmission of human t-cell lymphotropic virus types i and ii: association with viral load. JAIDS Journal of Acquired Immune Deficiency Syndromes, 12(2), 193–201. 10.1097/00042560-199606010-00014Search in Google Scholar PubMed

[25] Kawasuji, H., Takegoshi, Y., Kaneda, M., Ueno, A., Miyajima, Y., Kawago, K., …, Yamamoto, Y. (2020). Transmissibility of COVID-19 depends on the viral load around onset in adult and symptomatic patients. PloS One, 15(12), e0243597. 10.1371/journal.pone.0243597Search in Google Scholar PubMed PubMed Central

[26] Khan, M. A., Ali, Z., Dennis, L., Khan, I., Islam, S., Ullah, M., & Gul, T. (2015). Stability analysis of an SVIR epidemic model with non-linear saturated incidence rate. Applied Mathematical Sciences, 9(23), 1145–1158. 10.12988/ams.2015.41164Search in Google Scholar

[27] Leontitsis, A., Senok, A., Alsheikh-Ali, A., Al Nasser, Y., Loney, T., & Alshamsi, A. (2021). Seahir: A specialized compartmental model for COVID-19. International Journal of Environmental Research and Public Health, 18(5), 2667. 10.3390/ijerph18052667Search in Google Scholar PubMed PubMed Central

[28] Li, C., Xu, J., Liu, J., & Zhou, Y. (2020). The within-host viral kinetics of sars-cov-2. bioRxiv. 10.1101/2020.02.29.965418Search in Google Scholar

[29] Mandale, R., Kumar, A., Vamsi, D., & Srivastava, P. K. (2021). Dynamics of an infectious disease in the presence of saturated medical treatment of Holling type iii and self-protection. Journal of Biological Systems, 29(2), 245–289. 10.1142/S0218339021400064Search in Google Scholar

[30] Martcheva, M., & Li, X.-Z. (2013). Linking immunological and epidemiological dynamics of HIV: the case of super-infection. Journal of Biological Dynamics, 7(1), 161–182. 10.1080/17513758.2013.820358Search in Google Scholar PubMed PubMed Central

[31] Mellors, J. W., Rinaldo, C. R., Gupta, P., White, R. M., Todd, J. A., & Kingsley, L. A. (1996). Prognosis in hiv-1 infection predicted by the quantity of virus in plasma. Science, 272(5265), 1167–1170. 10.1126/science.272.5265.1167Search in Google Scholar PubMed

[32] Miller, E., Warburg, A., Novikov, I., Hailu, A., Volf, P., Seblova, V., & Huppert, A. (2014). Quantifying the contribution of hosts with different parasite concentrations to the transmission of visceral leishmaniasis in Ethiopia. PLoS Neglected Tropical Diseases, 8(10), e3288. 10.1371/journal.pntd.0003288Search in Google Scholar PubMed PubMed Central

[33] Murillo, L. N., Murillo, M. S., & Perelson, A. S. (2013). Towards multiscale modeling of influenza infection. Journal of Theoretical Biology, 332, 267–290. 10.1016/j.jtbi.2013.03.024Search in Google Scholar PubMed PubMed Central

[34] Mwalili, S., Kimathi, M., Ojiambo, V., Gathungu, D., & Mbogo, R. (2020). Seir model for COVID-19 dynamics incorporating the environment and social distancing. BMC Research Notes, 13(1), 1–5. 10.1186/s13104-020-05192-1Search in Google Scholar PubMed PubMed Central

[35] Ndaïrou, F., Area, I., Nieto, J. J., & Torres, D. F. (2020). Mathematical modeling of COVID-19 transmission dynamics with a case study of Wuhan. Chaos, Solitons & Fractals, 135, 109846. 10.1016/j.chaos.2020.109846Search in Google Scholar PubMed PubMed Central

[36] Prakash, D. B., Vamsi, D., Rajesh, D. B., & Sanjeevi, C. B. (2020). Control intervention strategies for within-host, between-host and their efficacy in the treatment, spread of COVID-19: A multi scale modeling approach. Computational and Mathematical Biophysics, 8(1), 198–210. 10.1515/cmb-2020-0111Search in Google Scholar

[37] Quinn, T. C., Wawer, M. J., Sewankambo, N., Serwadda, D., Li, C., Wabwire-Mangen, F., …, Gray, R. H. (2000). Viral load and heterosexual transmission of human immunodeficiency virus type 1. New England Journal of Medicine, 342(13), 921–929. 10.1056/NEJM200003303421303Search in Google Scholar PubMed

[38] Roldan, E. Q., Biasiotto, G., Magro, P., & Zanella, I. (2020). The possible mechanisms of action of 4-aminoquinolines (chloroquine/hydroxychloroquine) against sars-cov-2 infection (COVID-19), A role for iron homeostasis? Pharmacological Research, page 104904. 10.1016/j.phrs.2020.104904Search in Google Scholar PubMed PubMed Central

[39] Samui, P., Mondal, J., & Khajanchi, S. (2020). A mathematical model for COVID-19 transmission dynamics with a case study of India. Chaos, Solitons & Fractals, 140, 110173. 10.1016/j.chaos.2020.110173Search in Google Scholar PubMed PubMed Central

[40] Sarkar, K., Khajanchi, S., & Nieto, J. J. (2020). Modeling and forecasting the COVID-19 pandemic in india. Chaos, Solitons & Fractals, 139, 110049. 10.1016/j.chaos.2020.110049Search in Google Scholar PubMed PubMed Central

[41] Tu, Y.-F., Chien, C.-S., Yarmishyn, A. A., Lin, Y.-Y., Luo, Y.-H., Lin, Y.-T., …, Chiou, S.-H. (2020). A review of SARS-CoV-2 and the ongoing clinical trials. International Journal of Molecular Sciences, 21(7), 2657. 10.3390/ijms21072657Search in Google Scholar PubMed PubMed Central

[42] Van den Driessche, P. (2017). Reproduction numbers of infectious disease models. Infectious Disease Modelling, 2(3), 288–303. 10.1016/j.idm.2017.06.002Search in Google Scholar PubMed PubMed Central

[43] Wang, B. X., & Fish, E. N. (2019). Global virus outbreaks: Interferons as 1st responders. In: Seminars in Immunology (vol. 43, p. 101300). San Diego: Academic Press.10.1016/j.smim.2019.101300Search in Google Scholar PubMed PubMed Central

[44] Wang, T., Wu, Y., Lau, J. Y.-N., Yu, Y., Liu, L., Li, J., …, Jiang, B. (2020). A four-compartment model for the COVID-19 infection-implications on infection kinetics, control measures, and lockdown exit strategies. Precision Clinical Medicine, 3(2), 104–112. 10.1093/pcmedi/pbaa018Search in Google Scholar PubMed PubMed Central

[45] Wang, X., Wang, S., Wang, J., & Rong, L. (2022). A multiscale model of COVID-19 dynamics. Bulletin of Mathematical Biology, 84(9), 99. 10.1007/s11538-022-01058-8Search in Google Scholar PubMed PubMed Central

[46] Wen, W.-H., Chang, M.-H., Zhao, L.-L., Ni, Y.-H., Hsu, H.-Y., Wu, J.-F., …, Chen, H.-L. (2013). Mother-to-infant transmission of hepatitis b virus infection: significance of maternal viral load and strategies for intervention. Journal of Hepatology, 59(1), 24–30. 10.1016/j.jhep.2013.02.015Search in Google Scholar PubMed

[47] Zeb, A., Alzahrani, E., Erturk, V. S., & Zaman, G. (2020). Mathematical model for coronavirus disease 2019 (COVID-19) containing isolation class. BioMed Research International, 2020, 3452402. 10.1155/2020/3452402Search in Google Scholar PubMed PubMed Central

[48] Zhao, Z.-Y., Zhu, Y.-Z., Xu, J.-W., Hu, S.-X., Hu, Q.-Q., Lei, Z., …, Luo, L. (2020). A five-compartment model of age-specific transmissibility of SARS-CoV-2. Infectious Diseases of Poverty, 9(1), 1–15. 10.1186/s40249-020-00735-xSearch in Google Scholar PubMed PubMed Central

Received: 2023-12-21

Revised: 2024-05-06

Accepted: 2024-05-06

Published Online: 2024-06-21

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/cmb-2024-0006

Keywords for this article

multi-scale model; well-posedness; elasticity; comparative effectiveness; heat plot

Creative Commons

BY 4.0