The ‘wicked problem’ of telerehabilitation: Considerations for planning the way forward

1.   Motivations and plan of the paper

The emergence of the SARS-CoV-2 virus, responsible for the initial COVID-19 outbreak and the subsequent pandemic ^[3], has affected our minds, health, and well-being to the point of visibly altering the way we think and organize our lives, including our individual feelings, work organization, economics, and all expressions of collective behavior in our societies. These enormous problems have affected all countries of the world.

On the positive side, however, there has been a growing awareness that science is a primary asset to be respected and preserved. In addition, human societies are now aware that we live in a complex and interconnected world. Accordingly, science has generated an enormous amount of research activity, not only in the various branches of biology and medicine, such as virology, immunology, and medical care, but also in parallel fields such as economics ^[23,36], psychology, and social sciences ^[29]. Thus, a new vision of social organization and a greater attention to welfare problems is developing in our societies. A dialog within the scientific community was quickly opened, see ^[41], and brought to the attention of society the fragility of our planet ^[6].

The mathematical sciences have followed this activity, as can be seen in ^[8] and in the contents of various special issues on this topic ^[9,68]. Most of the published studies focus on various modifications of the susceptible-infectious-recovered (SIR) model, to incorporate specific aspects that could impact viral transmission; see ^[16,30,55]. Some applications refer specifically to regional areas ^[14,28]. The structural simplicity of these models is discussed, for example, in ^[21]. An important problem addressed by compartmental models of SIR-type is the calibration of the models with empirical data focusing on regional or national areas. These studies have to deal with the difficulty of collecting reliable empirical data, somewhat related to contagion problems related to transport dynamics.

On the other hand, some innovative ideas have been proposed by various authors to account for the different types of heterogeneity of human societies and their aggregations. For example, in ^[8] the authors proposed a multiscale vision related to a systems approach that goes far beyond deterministic population dynamics. Indeed, individual responses to infection and pandemic events are heterogeneously distributed throughout the population. In addition to biological heterogeneity, there are other types of heterogeneity, such as physical, related to age, and social heterogeneity that can include also the level of education, which leads to a greater or lesser ability to reach an advanced level of awareness regarding the risk of contagion ^[2,40]. In particular, contagion awareness also plays an important role in open areas and has led to various studies of crowd dynamics with internal behavioral variables ^{[1,12,13,18,39]}.

Another key issue in regard to SARS-CoV-2 infections is the in-host dynamics following viral infection, which can end with immediate recovery, hospitalisation or even death of patients. Experimental research over the past four years has led to a large amount of data on viral loads or various immune responses (see Figure 1 for an example of two such datasets). However, due to the heterogeneity in the data (e.g., viral RNA vs. sgRNA vs. virus titers ^[31]) and in the experiments, the difficulties in interpreting such data in the context of single-scale ODE models that consider some generic "viral load" counted from the moment the virus enters the body (and not from when symptoms appear, as shown in Figure 1 below), or even due to the complexity of the multiscale mathematical models that requires new mathematical techniques to parametrise the models across multiple scales (see the Discussion section below), not all of this data is currently being used to parametrise the mathematical models for in-host viral dynamics.

Figure 1.  Example of viral load data (viral copies/ml) and immune data (IFN$ \gamma $-secreting T cells, as measured by ELISPOT assays) taken from Challenger et al. ^[19]. The authors in ^[19] have summarized longitudinal viral load data from 19 different publications. They have also collected detailed data on immune responses from the study by Tan et al. ^[67]. Here we show this data for T cell numbers. Note the time axis showing "days from symptoms onset", which is different from the approach taken by the mathematical models that usually consider time from the start of infection.

DownLoad: Full-Size Img PowerPoint

There are currently a few publications that focus on modeling the multiscale heterogeneous aspects of in-host infections with SARS-CoV-2 ^[34,74]. For example, we mention the multiscale modeling study in ^[74], which considered the heterogeneous cells with different numbers of surface receptors (thus at the microscale level). The multiscale modeling study in ^[34] used an agent-based model to describe the heterogeneous cell population (where cells have different internal concentrations of an antiviral drug) and the heterogeneous cytokine concentrations secreted by the cells. These published models focused on discrete aspects of the variable that gives the heterogeneity of a population, and thus they were mainly described by ordinary differential equations.

This paper is devoted to modeling the multiscale in-host dynamics that develop within each individual immediately after viral infection. The multiscale aspect is related to the heterogeneity of the immune response to the viral infection (i.e., to viruses and infected epithelial cells), which determines the outcome of the infection. In contrast to previous studies that focus on discrete states of cells, here we assume a continuous activity variable that describes the immune response. This assumption is based on experimental studies that show a continuum of immune cell activation states, that results from various cell integration of cytokine signals that drive their activation ^[25,32]. Throughout this theoretical study, we consider the competition between the virus and the heterogeneous immune responses that appear immediately after the viral infection of the upper respiratory tract, before the virus moves to the lungs and causes further tissue damage.

The paper is structured as follows:

● Section 2 presents the biological assumptions regarding the dynamics of the virus versus the immune system, where the initial competition involves the innate immunity, while the adaptive immunity acts after a lapse of time necessary to learn the presence of the pathogen. The key feature of the dynamics is to quantify the outcome of the competition, i.e., to understand whether the innate immunity succeeds in preventing the virus from reaching the lungs, while the activation of the adaptive immunity by vaccines can avoid this event.

● Section 3 presents the general multiscale modeling framework that uses the so-called kinetic theory of active particles ^[10]. Then, this general framework is applied to a specific case of in-host viral infections of epithelial cells, which trigger heterogeneous immune anti-viral responses, where immune cells evolve different levels of activation. The details of the model (and the biological mechanisms behind different virus-cell and cell-cell interactions) are also specified in this section.

● Section 4 presents some numerical simulations that illustrate the range of model dynamics: from virus elimination, to virus re-infection, to virus persistence and immune blow-up. Different parameters are varied to explore the possible outcomes of the model: the mechanisms behind short-term or long-term infections, the impact of different initial viral loads, and the impact of persistant immunity.

● Finally, Section 5 looks ahead to new research perspectives by focusing on the potential of this multiscale modeling framework to understand the complex dynamics of heterogeneous immune responses.

2.   Phenomenology of the in-host dynamics

This section presents a phenomenological interpretation of the specific biological features of the in-host competition between the virus and the immune system ^[50], where the dynamics take place after contagion. As will be shown in the next section, this framework is intended to contribute to the derivation of mathematical models that can describe the aforementioned competition in terms of a differential system. We refer the reader to ^[8,44,47,62] for more details regarding the description of the in-host virus-immune competition.

Here, we focus on the early stages of infection, when the virus succeeds in bypassing the mucosal barrier of the upper respiratory tract and begins to infect the epithelial cells, where it is confronted by various immune reaction mechanisms. While innate responses play a very important role in blocking the initial stages of invasion and in initiating antiviral immunity, the adaptive immune response stops viral replication, leading to recovery from COVID-19 and the induction of persistent immune memory ^[44]. Unfortunately, sometimes the immune mechanisms fail to contain the invasion and the viral particles are able to reach the lung tissue, a critical factor in determining the evolution of the disease and the fate of the patient ^[47]. This can happen because the number of infectious virus particles is so high that they overcome the resistance of the immune mechanisms ^[73]. Indeed, the cumulative number of virus particles infecting a patient is not easy to estimate. It could range from a minimum number of virus particles needed to establish infection (i.e., infectious dose) to higher numbers that can be acquired in hospitals, crowded environments, or from super-spreading events. Alternatively, virus spread could be promoted by individual genetic variants ^[4], immunodeficiencies, even those that would otherwise have gone unnoticed ^[48], senescence ^[42], or comorbidities that impair the function of the immune system.

Referring to the study in ^[8], we remark that the investigation of viral infections in general requires a multiscale approach. The macroscale corresponds to individuals that may be infected or uninfected, while the microscale corresponds to in-host entities within infected individuals. The link between the two scales is provided by the dynamics of contagion, which depend on the multiple parameters mentioned above. Keeping this general framework in mind, the content of our paper focuses on the in-host dynamics evolving from the viral load as an input data. Also related to input data, we remark that since SARS-CoV-2 is a new virus ^[47], the antibodies induced by other coronaviruses, such as SARS-CoV and MERS, in most of cases do not recognize and neutralize this new virus.

The specific features to be considered in the modeling approach of the dynamics of the biological system under consideration are described by the following points, which are preliminary to modeling. We do not naively claim that the description is exhaustive, as it is limited to the specific assumptions that will be actually included in the modeling approach.

1. Modeling Framework: The modeling approach is developed at the micro-scale, within an infected individual, and focuses on the dynamics of the in-host entities, i.e., cells within the tissue, virus particles, and immune agents. After a viral infection, there is a competition between the virus and the immune system within each infected individual. The infection may progress (or regress) due to a superiority (or inferiority) of the virus over the immune defense, and may end with a full recovery of the patient, with the so-called long COVID-19, or with the death of the patient. These outcomes are based on a complex dynamic by which antiviral immune mechanisms prevent or not the passage of the virus from the upper respiratory tract to the lungs, as well as on the time needed for the activation of adaptive immunity to eliminate the virus.

2. Multiple Characteristics: As noted above, the macroscopic scale corresponds to individuals and the microscopic scale corresponds to in-host (i.e., within an individual) entities such as epithelial cells, virus particles, and immune agents.

3. Contagion: The likelihood of infection depends on the number of viral copies (viral load) that can reach the respiratory tract of a susceptible person. The arrival of an infectious viral load depends mostly on the dynamics at the macroscopic scale, i.e., the interaction between a person spreading the virus and other susceptible individuals and the social distance between them. The risk of contagion is also influenced by the characteristics of the area where the infection occurs (temperature, humidity, etc.) and on the time during which susceptible individuals are exposed to the virus. The result of these dynamics modulates the number of infectious virus copies, a modulation that critically influences the initial in-host dynamics of the virus-immune system within each susceptible individual.

4. Viral Load: Viral load is the number of virus particles simultaneously penetrating a person. The total number of SARS-CoV-2 virions was estimated in a recent experimental study ^[64]. The minimum infectious dose of SARS-CoV-2 that can cause COVID-19 in humans is still an open question, although some modeling studies and some experimental animal studies have estimated the minimum infectious dose of SARS-CoV-2 to be of the order of hundreds of virions ^[37,69]. In a recent experimental study, the total number of SARS-CoV-2 present in an infected individual has been estimated to be in the billions ^[64]. The highest number of virus particles is produced just before the onset of symptoms ^[37].

5. Affinity: Another important factor to consider is the affinity of the virus for human epithelial cells. Affinity is a specific property of each variant of the virus.

6. Dynamics of the Virus: Through spike proteins, the virus interacts with and enters human epithelial cells, where thousands of new copies of the virus are produced. These then invade additional epithelial cells.

7. Virus Variants–Copying Errors: During the rapid production of new viruses within the infected cell, many mistakes are made in the copying (transcription) of the viral RNA. Most of these errors (mutations) are inconsequential or result in the production of less effective virus particles. In a few cases, new variants with more invasive properties may appear. Although the likelihood of the emergence of more infectious variants is very low, Darwinian competition among virus particles may favor their spread.

8. Dynamics of the Immune System: On the cell membrane as well as in the cytoplasm, human cells have several sensors (e.g., pattern recognition receptors) that detect the presence of viral RNA. These sensors trigger a variety of reactions that lead to the destruction of viral RNA and the release of alarm signals. Among the early signals released by infected cells, a large family of interferons triggers cell antiviral activities and prevents the spread of viral infection to neighboring cells ^[52]. The virus-infected cell assembles the inflammasome and releases multiple pro-inflammatory cytokines. This triggers an initial localized inflammatory response. The infected cell may also undergo apoptosis before becoming a factory of new viral particles. On the other hand, several proteins encoded by the viruses can sabotage and dampen immune cell responses. The dynamics of the immune system's reaction to viral invasion are markedly affected by individual genetic variants that regulate the speed and intensity of each of the above-mentioned reaction mechanisms.

9. Activation of Adaptive Immunity: In most cases, adaptive immunity is not activated during these early stages unless primed by similar previous infections or by vaccines. However, at a slower pace, the presence of viral particles and viral proteins is perceived by adaptive immunity lymphocytes, alerted by the alarm signals released by both the infected epithelial cells and the innate immunity cells recruited into the initial local inflammatory response. Antibodies of different types (classes) are produced, progressively increasing in number and precision (affinity) toward the viral proteins.

The above biological model provides a basic scheme on which to derive the mathematical model. Some additional details will have to be added when a more precise definition of the biological events becomes important, as it will modulate the interpretation of the virus-immune system relationship. In Figure 2, the sequential steps of this biological model are shown: the dark red arrows show the progression of viral invasion, while the circles show the virus (V 1–4) and the immune defense activities (D1–7) that interfere with this progression. These interferences can be so important as to block infection and lead to viral clearance. Conversely, such perturbations can make viral progression unstoppable and catastrophic. The specific features of these events can be formalized by deriving a mathematical model.

Figure 2.  Representation of the in-host dynamics of the virus and of the immune response.

DownLoad: Full-Size Img PowerPoint

3.   Derivation of in-host mathematical models

This section is devoted to the derivation of mathematical models responsible for representing the in-host dynamics by transferring the biological description presented in Section 2 into a system of differential equations consistent with the mathematical kinetic theory of active particles (KTAP). The model derivation considers the following three steps:

1. Subdivision of the system into populations, called functional subsystems, shortly FSs, which play an effective role in the game. The state of each FS is defined by a distribution function over the microscale state, called activity, which is expressed by the individual entities, called active particles, shortly a-particles, that belong to each FS.

2. Derivation of a general mathematical structure suitable to define the dynamics of the distribution function which acts as the dependent variable of the dynamical system. As it is shown in ^[10], this structure should also capture, as far as it is possible, the main features of living systems. The derivation of the mathematical structure should also include the dynamics of learning (i.e., cells can learn from their interactions with other cells/particles and can change their activity state).

3. Modeling interactions involving active particles and inserting them into the given mathematical structure to obtain specific mathematical models. According to ^[10], these interactions can be nonlocal, nonlinearly additive, and can generate proliferative and destructive events.

These steps will be treated in more detail in the next subsections, where a multiscale mathematical model for the interactions between infected epithelial cells, virus particles, and immune cells with different activity levels will be proposed and further investigated numerically.

3.1.   Functional subsystems, activity, and representation

The state of each FS is given by a time-dependent distribution function over the activity variable $ u $, for instance $ f_i = f_i(t, u) $, $ i = 1, ..., \bar{N}_{f} $ (where $ \bar{N}_{f} $ is the number of FSs). If the activity is not modified by interactions, then the activity has a constant value $ u \cong u_0 $ and the system behaves as a deterministic population. In this case, we use the notation $ n_i = n_i(t; u_0) $ to distinguish the deterministic FSs from the stochastic FSs.

The following notation uses lowercase Roman letters for activities that can be modified by interactions and Greek letters for parameters that may be heterogeneously distributed but are not modified by interactions. The modeling approach also takes into account two reference quantities, the critical number of virus particles that can occupy the upper respiratory tract before entering the lungs ($ N_{M} $) and the maximum number of host cells in the upper respiratory tract ($ N_{H} $) that can be infected, which are used to normalize certain parameters involved in the dynamics to ensure appropriate units; see Section 4.

Let us now consider, with reference to Figure 3, the deterministic and stochastic FSs that play the game. A simple approach is to identify the following FSs and the activities expressed by each of them:

Figure 3.  Caricature description of the variables that form the functional subsystems considered in this study. The virus particles and the infected epithelial cells are described by densities, while the immune cells are described by a distribution over the activity variable.

DownLoad: Full-Size Img PowerPoint

$ i = 1 $: Viral particles express the ability to infect epithelial cells, where replication occurs. Since there is no progression of virus particles through an "activity-variable", the virus population can be considered deterministic, and the state of this 1-FS can be described by the density $ n_1 = n_1(t) $.

$ i = 2 $: Infected epithelial cells are the result of the viral infection spreading through the epithelial cell population. Here, we focus only on the infected cells (in the upper respiratory tract) which are less than the maximum number of host cells in the upper respiratory tract ($ N_{H} $) mentioned above. Again, since there is no progression of infected epithelial cells through an "activity-variable space", we can consider this cell population as deterministic and thus, the density of infected cells in this 2-FS is given by $ n_{2} = n_{2}(t) $.

$ i = 3 $: Immune cells learn the presence of the virus and try to control the infection via anti-viral and pro-inflammatory molecules, such as the interferons (IFN) ^[59]. For simplicity, here we average both innate and adaptive immune responses into a single variable that describes a general antiviral immune response. The state of this functional subsystem (3-FS) is given by the distribution $ f_3 = f_3(t; u). $ Here, $ u\geq0 $ describes the activation level of these immune cells (which is the result of different types and intensities of antigen exposure ^[33,51]), with $ u = 0 $ denoting the non-activated state, and $ u = 1 $ describing the fully activated cells.

Remark 3.1 Note that the variability in the activation level of immune cells (involved in both innate and adaptive immunity) is usually the result of variability in antigen exposure, with persistent antigen exposure leading to reduced immune cell function ^[33,51,72]. Therefore, as an example, this activity variable $ u $ could describe the state of exhaustion of CD8$ ^{+} $ T cells, as characterized by higher expression values for different surface molecules; e.g., high PD-1 expression (PD-1$ ^{HI} $), increased TCF-1 expression (from intermediate TCF-1$ ^{INT} $ values to high TCF-1$ ^{HI} $ values), or Granzyme B values (from low Granzyme B$ ^{LO} $ values to high Granzyme B$ ^{HI} $ values), as T cell exhaustion progresses ^[72]. This activity variable $ u $ could also describe the level of macrophage phagocytosis capability, as expressed by the surface FC receptor (as a loss in this receptor causes a reduction in the phagocytosis level of these innate immune cells) ^[56].

3.2.   Derivation of a kinetic theory mathematical structure

The derivation of mathematical models is based, as mentioned, on the kinetic theory of active particles. This statement can be made more precise by deriving a mathematical structure suitable to provide the framework to describe specific models. We consider the case of spatial homogeneity of models describing the time dynamics of the distribution functions:

where $ u $ is the microscale activity variable. The mathematical structure can be applied to all distribution functions, although some of the FSs correspond to a deterministic population. In fact, this case can be seen as a particular case of the most general one, where all dependent variables are stochastic.

The structure must include both conservative interactions, which modify the activity but not the number of active particles, and proliferative/destructive events, which also modify the number of particles. Guidelines for deriving such a structure are given in ^[10]. Accordingly, the said structure is as follows:

where $ f = \{f_1, f_2, f_3 \} $, where each component of $ f $ corresponds to a specific population. Accordingly, $ f_1 $ corresponds to viral particles (IP), $ f_2 $ corresponds to infected cells (FI), and $ f_{3} $ to immune cells (CI). The proliferative/destructive interaction functionals $ G_i [f] $ and $ L_i [f] $ could be local or non-local, constant or density-dependent, as we will see for a specific case in the next section. The progression probability $ \mathcal{A} $ could have different forms, as long as it is zero when $ u \approx 1 $ and it is one when $ u, u_{*}\approx 0 $; see also Remark 4.1. The interaction rate $ \eta_{ij} $ describes the frequency of interactions between population $ i $ and population $ j $ ^[10].

Remark 3.2 In the case of deterministic populations one has $ f_i(t; u_0) = n_i(t) \delta(u - u_0) $, where $ \delta $ denotes the delta distribution and $ n_{i}(t) $ denotes the density of the deterministic population $ i $. The structure (3.2) can be rapidly modified to include one or more deterministic populations, as shown by the model derived in the next subsection.

Remark 3.3 The density of deterministic population $ i $ is obtained by the zero-order moments of the distribution function as follows:

Remark 3.4 Since the innate immunity goes into action within minutes to hours of pathogen infection ^[45], we incorporate the effect of this innate immunity into the initial conditions of the model. Thus, $ f_{3}(t, u) $ will describe mainly the effect of adaptive immunity (with different activity levels $ u $). However, we need to be aware of this simplifying assumption (given that there is actually an overlap between the effects of innate and adaptive immunity).

3.3.   From interactions to derivation of models

The derivation of the mathematical model follows the assumptions on the phenomenological interpretation, in Section 2, of the biological system and of the description of the functional subsystems in Subsection 3.1. The variables and parameters that appear in the equations below are dimensional: variables are described in terms of cell densities and cell distributions (i.e., cell numbers/densities distributed over a range taken by the activity variable), as well as the number of viral RNA copies (RNAs$ /ml $ ^[7]), while parameters are described in terms of rates (e.g., proliferation rates associated with doubling times, and death rates associated with half lives) and carrying capacities (e.g., the maximum number of cells in the upper respiratory tract). For details on these parameters and variables, see Tables 1 and 2, in Section 4.

Table 1.  Initial conditions ($t = 0$) for numerical simulations of model (3.10).

Description	Initial value	Simulation range
Viral load	$n_{1}(t=0)=N_{0}= 10^{3}$	$[10^{2};10^{6}]$
Infected cells	$n_{2}(t=0)=1$	$[1; 10^{2}]$
Immune cells	$f_{3}(t=0, u)\;=10^{1}e^{-100u}$	$[e^{-100u}; 10^2 e^{-100u}]$

---

 | Show Table

DownLoad: CSV

Table 2.  Summary of parameters that appear in model (3.10), together with their dimensional units ("cells", "days", "RNA/ml"). The brackets in the 3rd column show the range over which the parameters are varied during the simulations. For the estimation of these values (and their scaling to ensure appropriate units), see the discussion in the main text.

Par.	Units	Baseline value; [range]	Description
$\alpha$	$\frac{ml}{day \times RNA}$	$1.44\times 10^{-6}$; $1.44\times [10^{-9}, 10^{-4}]$	rate at which virus infects epithelial cells [7]
$\gamma$	$\frac{1}{day}$	$0.69315$; $\left[0.1, 1\right]$	death rate of infected epithelial cells, associated with the production rate of new viral particles by these cells [54]
$\mu$	$\frac{1}{day\times cell}$	$10^{-6}$; $[10^{-4}, 10^{-10}]$	killing rate of virus particles and infected cells by anti-viral immune cells [27].
$\kappa(u)$	$\frac{ml}{day \times RNA}$	$10^{-9}$+$u 10^{-9}\;[10^{-11}, 10^{-6}]$	proliferation rate of immune cells in the presence of virus particles; we assume that cells with different activity levels $u\in [0, 1]$ have different proliferation rates [20,67]
$\lambda$	$\frac{1}{day}$	0.01 [$10^{-3}, 1$]	natural death rate of immune cells [17,20,53]
$\beta $	$\frac{1}{day}$	0.01 [$10^{-3}, 1$]	impact rate of virus particles on the activation status of immune cells
$N_{c}$	$\frac{RNA}{ml\times cell}$	$10^3\;[10^2, 10^4]$	the viral burst size, i.e., the number of virus particles released by an infected cell [7]
$\nu$	$-$	0.1 [$0, 1$]	coefficient describing the change in kernel $K(u)$ with respect to the activity variable $u$ (see Eq (3.6))
$N_{M} $	$\frac{RNA}{ml}$	$10^{9}$ [$10^{6}, 10^{11}$]	critical number of virus particles that can occupy the upper respiratory tract before invading the lungs
$N_{H}$	$cell$	$10^{5}$ [$10^{3}, 10^{9}$]	critical number of host cells in the upper respiratory tract that are infected [7]

---

 | Show Table

DownLoad: CSV

$ \bullet $ Dynamics of free viral particles ($ n_{1} $): Since the virus particles do not undergo conservative interactions to change their activity status, in Eq (3.2) we have $ \eta_{1j} = 0 $, $ \forall j\in\{1, 2, 3\} $. Regarding the proliferative interactions, we note that the viral particles are produced by the infected epithelial cells, and therefore in Eq (3.2) we have $ G_1 = \gamma N_{c}n_{2}(t) $. Here, $ \gamma $ is the death rate of infected cells that is associated with the release rate of new virus particles by these cells, and $ N_{c} $ is the number of virus particles inside one infected cell (see also Table 2). Regarding the destructive interactions, we note that the decay in the virion number is due to the action of immune cells which kill these viruses at a rate $ \mu $ (either directly through phagocytosis or indirectly through the secretion of cytokines such as IFN ^[59]). Thus, in Eq (3.2), we have:

where the kernel $ K(u) $ gives the activity action for the anti-viral effect of immune cells starting from the innate immunity corresponding to $ K = 1 $ to the increasing action of adaptative immunity. The following equation summarizes these biological assumptions:

Remark 3.5 $ K = K(u) $ is required to be an increasing function of $ u $. A simple model is as follows:

The heuristic model describes the growth of immune defence, where $ u = 0 $ denotes the the action of innate immunity only, while $ \nu\, u\uparrow \Rightarrow K \uparrow $ models the activation of the adaptative immunity.

Remark 3.6 Parameter $ \gamma $ is multiplied by the maximum viral load per cell $ N_{c} $ ^[7]. When infected cells burst open, all of these viral particles are released into the micro-environment (and contribute to the population of free virions). In ^[7], the authors estimated this burst size to be about $ 10^{3} $ virions, and thus we may take $ N_{c} = 10^3 $ virions/cell.

$ \bullet $ Dynamics of infected epithelial cells ($ n_{2} $): Since the progression of infected epithelial cells is not modified by interactions with other cells or viral particles, in Eq (3.2), we have $ \eta_{2j} = 0 $, $ \forall j\in \{1, 2, 3\} $. In regard to the proliferative interactions, we note that the number of infected cells increases, as in Eq (3.4) (at rate $ \alpha $) due to viral particles interacting with the uninfected cells ($ N_{H}-n_{2}(t) $) and entering these cells. The destructive interactions are the result of the anti-viral effect of immune cells ^[22,65], as well as due to virus-induced cell death at rate $ \gamma $ following viral proliferation (i.e., cells burst open and die). Thus, in Eq (3.2) we have

with $ K(u) $ a kernel defining the activity space for the anti-viral immune responses. These biological assumptions are summarized by the following equation:

where the abbreviation $ dvp $ denotes cell death due to virus proliferation.

$ \bullet $ Dynamics of immune cells ($ f_{3} $): As already mentioned, the distribution function $ f_3 $ includes both innate and adaptative immunity, where the activation of innate immunity corresponds to $ u = 0 $, while interactions with virus particles activate cells of adaptative immunity from $ u = 0 $ to $ u = 1 $. This dynamic is described by the probability density $ \mathcal{A}(u_{*}\to u) $ corresponding to the interaction rate $ \beta $ between immune cells and viral particles. In addition, these interactions induce a proliferation of immune cells (with rate $ \kappa $, which depends on the activation status $ u $). These assumptions lead to the following equation:

Therefore, the mathematical model consists in the following system of mixed ODEs and PDEs:

Remark 4.1 The progression probability $ \mathcal{A} $ could have different shapes, as long as it is approximately one when $ u, u_{*}\approx 0 $, it is zero when $ u = 1 $, and it is defined for $ u_{*}\geq u $. Examples are: (a) $ \frac{1-u}{1+(u-u_{*})^{+}} $ and (b) $ (1-u)e^{-a\frac{(u-u_{*})^{+}}{1-u}} $. For the simulations presented throughout this study, we use the function given in (a): $ \mathcal{A}(u_{*}\to u) = \frac{1-u}{1+(u-u_{*})^{+}} $.

4.   Numerical simulations

This section reports some sample simulations for model (3.10) with the purpose of exploring the qualitative behavior of the system considered here. We emphasise that while we tried to choose realistic parameter values for these simulations (see the discussion below on the initial conditions, and the discussion in Subsection 4.1 on all other model parameters), there is too much uncertainty in these values to even think about quantitative predictions and comparison with real observations (see Figure 1 and the discussion about model calibration at the end of Section 5).

The initial conditions used for these simulations are listed in Table 1. For consistency with the data in Figure 1, we assume that time represents the days from when symptoms appear. Since in ^[64] the authors estimated between $ 10^{2} $–$ 10^{4} $ virions in the nasal mucosa and pharynks, and the same numbers in the trachea and bronchus, as well as in the tonsils and lymph nodes of an infected individual, here we decided to start with $ n_{1}(t = 0) = 10^3 $ virions. In regard to the infected cells, in ^[64] the authors estimated between $ 10^{4} $–$ 10^{6} $ infected cells per person around the peak of the infection (under the assumption that 1 cell could be infected by 10 virions at a given time; more virions per cells would have been resulted in much lower estimations of infected cell numbers). These authors ^[64] also estimated around $ 10^{11} $ potential host cells, which resulted in $ 1 $ in $ 10^{5} $–$ 10^{7} $ cells that could be infected. Since we do not know how many cells can be infected in the upper respiratory tract (the focus of this study), we decided to start the simulations with $ n_{2}(t = 0) = 1 $ and varied these initial conditions between $ [1; 10^2] $. Finally, in regard to the immune responses, we note that the authors in ^[64] have estimated up to $ 10^9 $ SARS-CoV-2 specific T cells, with an unknown fraction found in infected tissues, and between $ 10^2 $–$ 10^4 $ T cells per infected cell. Since we do not know how many of these immune cells are fully activated (i.e., $ u $ close to 1), we assume an initial distribution $ f_{3}(t = 0, u) = 10e^{-100 u} $ and vary it between $ [10^{0}e^{-100 u}; 10^{2} e^{-100 u}] $; see also Table 1.

For the numerical simulations (implemented in C), we discretized the integrals in Eq (3.10) using Simpson's rule. For the progression probability, we used $ \mathcal{A}(u_{*}\to u) = \frac{1-u}{1+(u-u_{*})} $. Then, time was discretized using a Runge-Kutta method.

As the numerical simulations presented below will focus on the impact of various model parameters on the spread of viral infection through the epithelial cells, in the next sub-section, we discuss briefly the estimated values of these parameters (and the ranges over which they are varied).

4.1.   Estimating parameter values

The parameter values used for the numerical simulations in this study have been estimated as follows:

● In regard to the host cells (per person), in ^[7] the authors noted that there are approximately. $ 10^9 $ mucous cells in the nasal cavity. Throughout this theoretical study, we choose a baseline value $ N_{H} = 10^{5} $ much lower than the maximum number of host cells in the upper respiratory tract, since it makes sense to assume that not all cells will be infected before the virus invades the lungs. We also run simulations for $ N_{H}\in[10^{3}, 10^{9}] $.

● In regard to the critical number of virus particles that occupy the upper respiratory tract before invading the lungs, we note that the number of SARS-CoV-2 genomic copies in nasopharyngeal swabs is more than $ 10^{6} $ viral RNAs copies/ml when the symptoms appear ^[60], and persists for approximately 5 days before declining ^[60]. In ^[7], the authors noted that there are between $ 10^{6} $–$ 10^{11} $ RNAs$ /ml $ in sputum, $ 10^{6} $–$ 10^{9} $ RNAs/swab in the nasopharynx and $ 10^{4} $–$ 10^{8} $ RNAs/swab in the throat. Moreover, as mentioned above, in ^[64] the authors estimated between $ 10^{2} $–$ 10^{4} $ virions (or $ 10^{6} $–$ 10^{8} $ RNA copies) in the nasal mucosa and pharynx, and the same numbers ($ 10^{2} $–$ 10^{4} $ virions or $ 10^6 $–$ 10^9 $ RNA copies) in the trachea and bronchus, and again in the tonsils and lymph nodes of an infected reference individual with a total body weight of 70 kg. The same authors estimated between $ 10^5 $–$ 10^7 $ infectious units (virions) in lung tissue. Due to the variability in data across all of these different studies, here we choose a baseline value $ N_{M} $ = $ 10^{9} $, but we also run simulations for $ N_{M}\in[10^{6}, 10^{11}] $.

● In regard to the replication timescale of SARS-CoV-2, in ^[7] the authors noted that the time for a virion to enter a cell is approximately $ 10 $ min, which translates into an infection rate $ \approx 144/day $. However, since the units for $ \alpha $ are $ 1/(day \times RNAs/ml) $, we need to rescale the infection rate by virion levels. To this end, we choose $ \alpha\approx 144/N_{M} $. As a baseline value, we thus choose $ \alpha = 1.44\times 10^{-6} $, but we vary this parameter within the range $ \alpha \in [1.44\times 10^{-9}, 1.44\times 10^{-4}] $.

● In ^[54], the authors show experimentally that infected human cells reduce their viability by a half within 24 hours. From here, we estimate a death rate of $ \gamma = \ln(2)/1\: $day = $ 0.693/ $day.

● The burst size (i.e., number of virions released by one infected cell) can vary between $ 10 $–$ 10^{2} $ in ^[64] and $ 10^3 $ in ^[7]. Here, we choose $ N_{c} = 10^{3} $.

● Innate immune responses are observed within hours. Anti-viral adaptive T cell immune responses are observed within 7 days of symptoms ^[49]. In fact, the study in ^[67] showed the presence of T cell responses around days 3–5 from symptom onset in mild COVID-19 patients, and around day 10 in moderate/severe patients. Earlier quantitative studies on the kinetics of murine T cells in response to different viral infections ^[20] calculated the doubling time of T cells during the initial expansion phase to be between 8–11 h, which translates into a murine proliferation rate between $ 1.5 $–$ 2.0 $/day. In healthy humans, this proliferation rate is much slower ^[43] and was shown to differ between the naive T cells ($ 0.0005 $–$ 0.002/ $day), stem cell memory T cells ($ 0.0007 $–$ 0.007 $/day), effector-memory ($ 0.042 $/day), and central memory ($ 0.01 $/day) T cells. It is expected that human viral infections lead to a faster immune proliferation rate. We can describe all these biological observations (showing an increased proliferation rate for more active cells, i.e., the effector cells) with a simple linear proliferation function $ \kappa(u) = 1+u $. Due to the units for the immune proliferation rates (i.e., $ ml/(day\times RNAs) $; see Table 2), we rescale the above expression by $ N_{M} $ and choose average proliferation rates $ \kappa(u) = (10^{-9} +u\times 10^{-9})ml/(day\times RNAs) $, with $ u\in [0, 1] $. However, we vary $ \kappa \in [10^{-11}, 10^{-6}]\:ml/(day\times RNAs) $ to account also for the variability in the proliferation of T cells depending on the viruses used.

● In ^[53], the authors mentioned that following their expansion in response to viral particles, the T cell population then decays with a half-life of approximately 200 days, which translates into a death rate $ \lambda = 0.0035 $/day. However, this includes also the memory response that persists much longer after the death of effector cells. An earlier study ^[17] (not focused on SARS-CoV-2) estimated that effector memory T cells can have a life span from 6 days to 6 weeks. Another very early study ^[20] calculated the half-life of murine T cells during the contraction phase following the peak of the response between 41 h and 3 days, which translates into a cell death rate $ \lambda\in(0.23, 0.4) $. Since different viruses can trigger different immune responses with different cell kinetics (for either murine or human cells), throughout this study we consider a baseline value $ \lambda = 10^{-2} $, but we vary this rate within a range $ \lambda\in(10^{-3}, 10^{1}) $. This large parameter range is the result of the fact that in this study, we do not distinguish between the different types of immune cells (e.g., effector or effector-memory T cells, or the new concept of memory macrophages ^[70]) which have completely different lifespans described by different $ \lambda $.

● Since the COVID-19 patients usually recover within 10–14 days, we can assume that within 10 days the virus is eliminated by the immune cells. Due to the units of $ \mu $ (i.e., $ 1/(day\times cell) $), we rescale by $ N_{H} $ and thus we take an average baseline value $ \mu \approx 10^{-6} $. Note that experimental studies to quantify the efficiency of human CD8 T cells against specific viruses (thus also against SARS-CoV-2) are non-trivial ^[27], and previous estimates of this killing rate were performed on animal cells infected with different viruses, which led to killing rates between $ 0.3 $–$ 0.8 $/day ^[27]. Using all of this information, in this study, we consider the following range for the killing rate of viruses and infected cells by the immune cells: $ \mu\in(10^{-10}, 10^{-4}) $.

4.2.   Model dynamics

In the following, we illustrate numerically some of the dynamics of kinetic (multiscale) model (3.10).

● Baseline dynamics. We start our numerical investigation into the behavior of model (3.10) by showing in Figure 4 the baseline dynamics of this model. The simulations are performed with the baseline parameter values from Table 2 and the initial conditions from Table 1. We see that, for these parameter values, the virus particles and the infected epithelial cells are eliminated after the first two days following the onset of symptoms. The elimination is the result of an increase in the level of immune cells. Note that experimental/clinical studies have observed a rapid drop in viral load after the peak of the infection ^[19,64]; see also Figure 1. Moreover, in Figure 4 we notice a prolonged immune response, which also seems to be consistent (at least qualitatively) with clinical data shown in ^[19] and Figure 1.

Figure 4.  Baseline dynamics of continuum model (3.10): (a) density of virus particles: $ n_{1}(t) $; (b) density of infected cells: $ n_{2}(t) $; (c) distribution of immune cells across the activity space: $ f_{3}(t, u) $. All parameters are as in Table 2. In (a) and (b), the sub-panels (i) show the 2D plots $ n_{1, 2}(t) $ vs. time, while the sub-panels (ii) show the 3D plots $ n_{1, 2} $ vs. $ u $ vs. time. The inset (right) sub-figure in (c) shows the same graph as in the left main figure, but for time $ t\in[0, 0.8] $, so that we see more clearly how the immune cells start to become activated (i.e., larger $ f_{3}(u, t) $ as $ u\to 1 $ and $ t\to 0.8 $).

DownLoad: Full-Size Img PowerPoint

However, for the parameter values considered here, this viral elimination is not permanent and the infection can relapse in the very long term. Such a re-infection behavior is consistent with various clinical studies ^[35,63,66], which noted that the rates of re-infections were variable (from 3.1%–13.0% in ^[66], to 47% in ^[63]), and they could occur at different times (e.g., in ^[66] re-infection occurred approximately 180 days after the initial infection; in ^[63] re-infection occurred between 340–987 days after the initial infection). In Figure 5a, we show that for the baseline parameter values listed in Table 2, the infection returns after approximately 600 days. However, the second time there are less viral particles, less infected cells, and even a weaker immune response. This time period between the infection/re-infection depends on the persistance of functional immune cells that can recognise the virus and eliminate it (so model (3.10) incorporates implicitly an anti-viral memory immune response). In Figure 5b we investigate what happens when we reduce the persistence of functional immune cells (with different activation levels) by increasing their death rate $ \lambda $ to $ \lambda = 0.05 $. We see that in this case the second infection is slightly stronger than the first one (i.e., higher numbers of virus particles and infected cells, and higher levels of immune responses), while the third infection is slightly stronger than the second one. Moreover, for larger $ \lambda $ the time between the first and second infection is larger than the time between the second and third infection. We conclude the description of model dynamics shown in Figure 5 by noting that more than one re-infection is possible, as observed in some clinical studies ^[63].

Figure 5.  (a) Long-term baseline dynamics of model (3.10); the parameters are as in Figure 4. (b) Long-term dynamics obtained when we increase the death rate of immune cells to $ \lambda = 0.05 $; all other parameters are as in Table 2. In the sub-panels we show: (i) a 2D plot of the density of virus particles $ n_{1}(t) $; (ii) a 2D plot of the density of infected cells; and (iii) 3D plot of the distribution of immune cells $ f_{3}(u, t) $. We see that, for (a) the second infection is weaker than the first one, while for (b), the second infection is stronger than the first one.

DownLoad: Full-Size Img PowerPoint

● Increasing immune cell death. To have a better understanding of the impact of parameter $ \lambda $ (that incorporates implicitly memory immune responses to the virus), in Figure 6 we compare the characteristics of the model dynamics as we vary $ \lambda\in[0.01, 1] $. In panel (a), we see that the time between the original infection and the first re-infection decreases exponentially as we increase $ \lambda $: very fast for small $ \lambda $ (when there are enough immune cells in the system), and then much slower for large $ \lambda $ (when the immune cells are eliminated). Since, in Figure 5, we have seen that the second viral load could be higher or smaller than the first viral load, in panel (b) we show how the difference in the maximum viral load between the first and second infections varies with respect to $ \lambda $. The threshold between positive and negative values occurs at $ \lambda = 0.0275 $, and this corresponds to a delay between the 1st and 2nd infection of approximately 400 days (see the inset sub-figure in Figure 6a), which is consistent with clinical observations ^[63]. This result suggests that if the immune response is strong enough and persists for a very long time (i.e., when $ \lambda < 0.0275 $) then the second wave of re-infection will be smaller than the first wave. However, if the immune response is weak and does not last long (i.e., $ \lambda > 0.0275 $), the second wave of re-infection is higher than the first wave; but this difference between the first and second wave does not increase as $ \lambda $ increases significantly toward $ \lambda = 1 $ or past this value.

Figure 6.  (a) Time between the $ 1^{st} $ and $ 2^{nd} $ viral infections, as we vary $ \lambda $. (b) Difference in the maximum viral load during the $ 1^{st} $ and $ 2^{nd} $ viral infections: the threshold between positive values (i.e., $ 1^{st} $ wave higher than the $ 2^{nd} $ wave) and negative values (i.e., $ 1^{st} $ wave lower than the $ 2^{nd} $ wave) is at $ \lambda\approx 0.0275 $. The inset sub-figures are the same as the main figures, but with the horizontal axis on a logscale. All other model parameters are as in Table 2.

DownLoad: Full-Size Img PowerPoint

● No anti-viral immune response (i.e., immune cell exhaustion). We also investigate numerically what happens if there is almost no immune response against the virus or infected cells (and thus we take $ \mu = 0 $), while we keep the assumption that the immune cells can still be activated by the presence of the pathogen (i.e., $ \kappa(u) > 0 $) but they do not eliminate it. This assumption corresponds to clinical observations of high levels of immune cell activation and also high levels of immune cell exhaustion in hospitalized patients with very severe COVID-19 infections (some patients recovering, some dying) ^[71]. In Figure 7, we see that in this case, there is an over-activation of the immune response, which then leads to a blow-up in the solution (i.e., the blow-up is in the immune response). This blow-up is caused by the higher proliferation/activation rate of immune cells in the presence of the persistent pathogen (that also keeps proliferating).

Figure 7.  Model dynamics when we assume that the immune cells cannot eliminate the virus and the virus-infected epithelial cells: here $ \mu = 0 $. All other parameters are as in Table 2. The inset sub-figures in panel (c) show the immune response characterized by different activation states immediately after viral infection, i.e., $ t\in(0, 0.8) $, and for $ t\in (2.9, 3) $ just before the blow-up that occurs on the fourth day-not shown here.

DownLoad: Full-Size Img PowerPoint

● Varying virus initial conditiosn $ n_{1}(0) $, $ N_{H} $ and $ N_{M} $. Since experimental studies ^[26] suggested that initial viral load could predict the severity of illness, in Figure 8a, we show the time-evolution of $ n_{1}(t) $ as we increase its initial value to $ n_{1}(0) = 10^{4} $. We see that for the parameter values in Table 2, the virus is eliminated faster (due to a faster activation of the immune cells, in response to the higher initial viral load-not shown here).

Figure 8.  (a) Dynamics of $ n_{1}(t) $ when we increase the initial viral load: from $ n_{1}(0) = 10^3 $ to $ n_{1}(0) = 10^4 $. All other parameters are as in Table 2. (b) Dynamics of $ n_{1}(t) $ when we decrease $ N_{H} $ (from baseline $ N_{H} = 10^5 $) and increase $ N_{M} $ (from baseline $ N_{M} = 10^9 $). For a more thorough investigation of the effect of varying $ N_{H} $ and $ N_{M} $ on viral load, in panel (c) we show a colormap of the time when $ n_{1} $ reaches a maximum value, as we vary $ N_{M}\in[10^{8}, 10^{11}] $ and $ N_{H}\in[10^{3}, 10^{6}] $. The inset sub-figures in panel (c) show: (i) a change in the time when $ n_{1}(t) $ reaches its maximum value for $ N_{H} = 10^{4} $ and $ N_{M}\in [10^{8}, 10^{11}] $; (ii) the 3D plots of $ n_{1}(t) $ vs. time $ (t) $ vs. activity $ (u) $ for 2 specific cases: $ N_{M} = 10^{10} $ (top) where max($ n_{1} $) is reached at $ t\approx 5.5 $ but the virus is eliminated after $ t = 8 $; $ N_{M} = 10^{11} $ (bottom) where max($ n_{1} $) is reached at $ t\approx 5 $ but the virus is eliminated after $ t = 9.5 $.

DownLoad: Full-Size Img PowerPoint

These results depend significantly on two important parameters which could not be estimated from the available literature: $ N_{M} $ (the critical number of virus particles that can occupy the upper respiratory tract) and $ N_{H} $ (the critical number of host cells in the upper respiratory tract). In Figure 8, we show one example of the effect of these parameters on viral evolution: for $ N_{H} = 10^{4} $ and $ N_{M} = 10^{10} $, the amount of virus in the system is much lower (max($ n_{1}(t) $) = 3000, obtained at $ t = 5.5 $), but the infection lasts longer (until $ t = 8 $). For a more comprehensive investigation of the effect of these two parameters on viral persistence, in Figure 8c we show a colormap for the time when $ n_{1}(t) $ reaches its maximum (during the first infection) vs. $ N_{M}\in[10^{8}, 10^{11}] $ vs. $ N_{H}\in[10^{3}, 10^{6}] $. In the inset (i) we show the same time when $ n_{1} $ reaches its maximum, as we fix $ N_{H} = 10^{4} $ (as in sub-panel (b)) and we vary $ N_{M}\in [10^{8}, 10^{11}] $. First, we note that increasing $ N_{M} $ delays the time when $ n_{1}(t) $ reaches its maximum value. Second, we see that for $ N_{M} = 10^{10} $ max$ (n_{1}) $ is obtained for $ t = 5.5 $, while for $ N_{M} = 10^{11} $ max$ (n_{1}) $ is obtained for $ t = 5 $. To clarify this aspect, in sub-panels (ii), we show two 3D plots $ n_{1}(t) $ vs. $ t $ vs. immune activity $ u' $ for: $ N_{M} = 10^{10} $ (top) and $ N_{M} = 10^{11} $ (bottom). It is clear that for $ N_{M} = 10^{10} $, the $ n_{1} $ peak is higher and is reached later, but the infection is eliminated faster (due to the immune response–not shown here). In contract, for $ N_{M} = 10^{11} $ the $ n_{1} $ peak is lower and is reached sooner, but the infection continues to persist a bit longer.

5.   Summary and conclusion

This paper deals with the modeling and simulation of the in-host dynamics of the COVID-19 viral infections in the upper respiratory tract. The topic aims to shed new light on the complexity of heterogeneous immune-virus interactions, that could have an impact on patients. In fact, if the invasion of viral particles succeeds in moving to the lungs, a patient needs intensive care and the risk of death is greatly increased. Because of this motivation, the British Royal Society has decided to promote, right from the beginning of the pandemic, a voluntary team, coordinated by Mark Chaplain, to study this specific topic ^[61].

The plan of this paper was first to define the key features of the biological dynamics, and then to develop some multiscale models whose dynamics could be explored numerically. We began this study by proposing a general multiscale kinetic framework that could be used to investigate the infection of healthy cells by a viral population. This framework, based on the kinetic theory of active particles (KTAP) ^[8,9,10], can account for the heterogeneous nature of the antiviral immune responses with a continuous of activation states ^[32]. The resulting model is a mixture of equations for changes in the density of cells and viruses (which do not depend on the internal "activity" variable) and an equation for changes in the distribution of immune cells with different activity levels described by a continuous "activity" variable.

This KTAP theoretical framework (which can be applied to all viral infections and antiviral immune responses) was then applied in detail to the case of SARS-CoV-2 infections. To this end, we identified possible cell-cell and cell-virus interactions based on published experimental studies of COVID-19 antiviral responses. We also estimated the model parameters based on the published literature. Nevertheless, we are aware that the estimated values incorporated a lot of uncertainty and the results do not have any predictive power. Even so, Section 4 shows the possible theoretical outcomes of the model dynamics: from virus elimination in the short/medium time (Figure 4), to virus relapse in the long term (Figure 5), and even to a blow-up in the immune response due to an over-activation in response to virus accumulation (Figure 7). Some of these results were qualitatively similar to clinical data presented in Figure 1 (i.e., the high viral load observed when symptoms appear ^[19], the immune response that persists much longer after virus elimination ^[19,67]), or other clinical data showing re-infection of individuals already tested positive ^[35]. Other results, such as the blow-up of the immune response, only showed the potential dynamics of this model (as we could not find immune cell data for such cases. Most of the publications focused on cytokine levels in the context of tissue damage and patient mortality ^[38]). We also investigated the impact on model dynamics when we changed the initial levels of viral particles (Figure 8a), or we change different model parameters: $ \lambda $, $ N_{H} $ and $ N_{M} $ (Figure 8b, c). The simulations showed the possibility of having lower/higher secondary infections, as well as shorter/longer durations for these infections, all depending on the values of these model parameters (not easily estimated). Moreover, these simulation results of viral re-infections were consistent with the numerical results of other types of models that equally show re-infections: from discrete individual-based models ^[57] to continuous differential equation models ^[5]. Nevertheless, for a more clear understanding of the mathematical and biological mechanisms responsible for this oscillatory behavior, as well as for all other type of behavior, exhibited by this kinetic model (e.g., density blow-up), we need to focus on a more detailed analytical investigation of this model. However, such an investigation is not the purpose of this study, and will form the subject of a future investigation.

We emphasise again that this model does not have predictive power, since model dynamics depends very much on the unknown microscale-level parameters and functions, such as $ \mu $, and $ \kappa(u) $, as well as on the unknown macroscale-level parameters such as $ N_{H} $ and $ N_{M} $, which were approximated for the sake of simulations as their exact values are unknown (even if some lower/upper limits can be found in various studies—see Table 2). Moreover, the immune defences are heterogeneously distributed in the population (but we do not have any empirical data for such heterogeneous distribution of cell phenotypes); furthermore, their dynamics depend on the initial (unknown) viral load, which can have an important influence on the quantitative behavior, of the dynamics inside the various tissues of the body in the upper and lower respiratory tracts. Due to all of these unknowns, we emphasize that this model can only explore the qualitative behavior, of the system under consideration. Although some parameters have been chosen based on empirical data (see Table 2), and although some simulation results are qualitatively similar to some real data (see Figure 1), this study is not focused on model calibration.

Calibration of multiscale models with empirical data is not an easy task ^[24], since (i) biological processes at different scales are interconnected, but data (usually collected at a single scale as a result of one single experiment) might not account for such interconnectivity; (ii) it is not always very clear how to estimate the parameters that connect the different scales. In addition, some of model parameters (such as those triggering transitions in cells' phenotypes) are likely stochastic, and therefore they could be known probably only up to a distribution. For such cases, sensitivity and uncertainty analysis could be employed to better understand the impact of these parameters on model outcomes ^[58]. A complete global sensitivity analysis will be the subject of a future theoretical study focused on a comparison between kinetic and deterministic ODE models (as sensitivity and uncertainty analyses are now classical approches for ODE models, but not for kinetic models).

In summary, we believe that further studies should explore this data calibration problem for such kinetic multiscale models, as these models are sufficiently flexible to take into account all the aforementioned multiscale features of the biological systems. Additional studies can further develop the results of our paper: for instance, by introducing the dynamics of virus variants and vaccination programs into this multiscale approach (see ^[11,12,46]), while accounting also for their impacts on the immune responses. This would allow us to connect these in-host multiscale models with other between-host (either single-scale or multiscale) models for the spread of the disease through a human/animal population, in the presence or not of vaccination, as done in ^[15].

Author contributions

Nicola Bellomo, Raluca Eftimie, and Guido Forni : conceptualisation and writing (original draft, as well as reviewed draft and editing).

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

Nicola Bellomo acknowledges the support of the University of Granada, Project Modeling in Nature MNat from micro to macro, https://www.modelingnature.org.

Raluca Eftimie acknowledges funding from the MODCOV19 platform of the National Institute of Mathematical Sciences and their Interactions (CNRS).

Conflict of interest

The authors declare that there are no conflict of interest.