Multiscale models of Covid-19 with mutations and variants

1.   Aims and plan of the paper

The onset of SARS–CoV–2, responsible for the initial COVID-19 outbreak, has generated a subsequent pandemic all over the globe and societies. A multiscale modeling approach has been developed in [7] that involves many disciplines. In addition to applied mathematicians, also immunologists, economists and virologists contribute to the new approach. We refer to [26] which, in addition to a detailed description of the biology of the virus dynamics, provides useful indications on the world spread of a pandemic which has induced huge problems affecting not only health and well being in general, but also complex economical problems and has modified collective behaviors and social life. A concise quotation extracted from [7] summarizes the aforementioned concepts:

SARS–CoV–2 is mainly transmitted via respiratory droplets that an infected person expels. If the viral charge is high, the carrier is more infective. The large spike protein forms a sort of crown on the surface of the viral particles and acts as an anchor allowing the virus to bind to the Angiotensin-Converting Enzyme 2 (ACE2) receptors on the host cell. After binding, the host cell transmembrane proteases cut the Spike proteins, allowing the virus surface to approach the cell membrane, fuse with it and the viral RNA enter the cell. The virus hijacks the cell machinery and the cell dies releasing millions of new viruses thus generating a virus infection.

The interested reader is referred to [1,38,44,49] for the origin and mutations of the virus and to [46] for the RNA virus evolution. Modifications of the so-called SIR models [33], for instance SEIR models, have been proposed to consider the specificity of COVID-19 [31]. An interesting alternative is developed in [27] based on individual-based Markov models. Space propagation on networks of nodes have been developed in [14], where the dynamics within each node are modeled by SEIR structures, while transport dynamics are modeled with a kinetic theory approach. The general framework is the multiscale systems approach proposed in [7].

Various aspects of the immune competition are treated in [40,49], a topic somehow related to the search for vaccines [16,29,45]. Additional topics on the biology of the virus are reported in [17,21,36,41,47]. Various aspects of contagion dynamics have been studied in [30]. In addition, the impact on the economy of nations is treated in various research contributions, e.g. [3,5,23,24,32]. These citations do not claim to be exhaustive. Indeed, the research activity in the field keeps growing, motivated by the impact that the virus spread is exerting on our society and a huge flow of new results reaches daily the literature in the field.

A key indication from the literature in biology, is that infectivity is not a constant parameter, but it depends on the viral charge. In addition, it is necessary to account for the dynamics of the in-host immune competition which develops as the presence of virus particles releases danger signals which activate the reaction of the host's innate immunity. Coronaviruses are successful at suppressing various mechanisms, but not all of them, in the immune response.

Accordingly, the modeling approach should refer to general issues of the immune response [42], specifically, of the activation of the system from the innate to the activated immunity [18,20]. Pioneering contributions on the mathematical modeling of the immune competition are known in the literature, for instance [11] on the interaction of immune and tumor cells, while [22] specifically focuses on virus epidemics in the presence of mutations and selection.

Recently published articles have shown a rapidly growing interest towards a systems approach to social dynamics and behavioral economy, where the mathematical sciences are charged with capturing the complex features of these systems under the influence of individual and collective human behaviors. The objective of our paper specifically focuses on this topic by further developing the model [7] with the aim of depicting, firstly how relaxations of the confinement rules can generate sequential waves, and subsequently of describing how this specific action can modify the dynamics of the pandemics. These objectives can show how mathematics may provide to crisis managers simulations describing the possible scenarios generated by mathematical models and how external actions modify them thus contributing to planning strategy of the aforementioned actions.

We do not naively claim that mathematics can solve problems of epidemiology and virology as we are aware that the modeling approach cannot be developed by standalone mathematicians, an interdisciplinary approach is necessary as shown in [7]. Bearing all above in mind a detailed description of the contents of our paper is given in the following.

Section 2 presents the derivation of a multiscale mathematical model which describes the onset and diffusion of a COVID-19 epidemic as a natural development of the model proposed in [7]. The new model includes a detailed description of the confinement strategy on the epidemic dynamics and considers mutations into new variants. A multiscale approach is developed to model these dynamics which is followed by a in-host competition between virus and immune cells.

Section 3 presents first a set of simulations which aims at describing the scenarios of onset of the virus diffusion and of the subsequent waves related to locking and de-locking dynamics modeled by the parameters $\alpha_\ell$ and $\alpha_d$, respectively. Simulations are developed in absence of mutations to describe the dynamics of infected individuals and of the progression of the pathology which might lead to full recovery or to death.

Section 4 presents a second set of simulations which includes the dynamics of mutations and Darwinist selection. The objective consists in studying how this dynamic modifies the progression of the pathology with respect to that delivered in Section 3.

Section 5 develops a critical analysis on the use of simulations and related interpretation as a support to crisis managing and, out of this analysis, looks ahead to research perspectives.

2.   On a multiscale model of virus pandemics

A new class of mathematical models is derived in this section within the framework of the rationale proposed in [7]. Accordingly, we go beyond the approach of deterministic population dynamics by considering various aspects of the heterogeneity of individuals and of their reaction to the infection. The modeling is developed by a multiscale approach which refers both to dynamics of contagion and to the in-host competition between virus particles and the immune system. The two scales constantly interact, as the contagion at the macro scale depends of the viral load of each individual which in turn depends on the dynamics at the micro-scale. After contagion, the in-host dynamics include the immune competition within each individual.

We consider also a detailed description of confinement rules, as well as virus mutations and selection up to new variants. The contents are presented through a sequence of subsections focused on a phenomenological description of the biological system and its mathematical representation, derivation of a general mathematical structure suitable to capture the main features of the system, followed by derivation of models by implementing into said structure the mathematical description of multiscale interactions. The modeling is grounded on mathematical tools of the so-called kinetic theory of active particles [8] which provides mathematical structures suitable to model systems of interacting living entities. Applications of kinetic theory methods to living systems are critically analyzed in [2], where the main focus refers to dynamics far from equilibrium.

2.1.   Phenomenological description and representation

The following qualitative description of the biological dynamics is proposed to be transferred into mathematical models.

1. Interacting entities are modeled as active particles which are carriers of an internal state, called activity. These are divided into sub-populations called functional subsystems, in short FSs.

2. Two scales are used, namely the macro-scale corresponding to individuals who might be infected or not-infected and the micro-scale modeling the in-host dynamics within infected individuals. The micro-scale description includes different levels of the proliferative properties of the virus and of the defence ability of the immune system.

3. Contagion probability depends on the level of the infection, namely on the so-called viral charge, as well as on the social distance between individuals which might depend on time.

4. Within each infected individual, a competition occurs between the proliferative virus and the immune system. The level of infection can progress (or regress) in time due to a prevalence (or lack of prevalence) of the virus aggressiveness over the immune defence ending up with full recovery or death.

5. Mutations and selection up to the onset of new variants of the virus are modeled by post-Darwinist dynamics see [44] for general topics, while specialized issues are treated in [1,21,38].

We consider a population of $N_0$ individuals, which includes, at $t = 0$, a small fraction $n_0$ of infected individuals, namely a-particles. The time $t$ is defined in a bounded interval $[0,T]$. All dependent variables, which represent the state of the system, are referred (divided) to $N_0$. The representation of the system, according to our approach, is as follows:

$\bullet$ The population is subdivided into five sub-populations labeled by the subscripts $i = 1, \ldots, 5$. The abbreviation $i$-FS is used to denote the $i$-th population viewed as a functional subsystem. The subscript denotes the following FSs: $i = 1$, healthy; $i = 2$, infected individuals; $i = 3$, recovered from the infection; $i = 4$, dead; $i = 5$, infected by a new variant.

$\bullet$ The micro-state of the FSs includes two variables $u \in [0,1]$ and $w \in [0,1]$ corresponding to the progression of virus invasion and to the level of activation of the immune defence, respectively. In more details, $u = 0$ represents the absence of the viral infection, $u>0$ characterizes the presence of the disease; while $w = 0$ and $w = 1$ correspond, respectively, to the lowest (corresponding to the sentinel level), and highest activation immune system.

$\bullet$ Increasing values of $u$ and $w$ correspond to increasing values of proliferative activities of the virus and of the immune system. Discrete variables are used:

where the same number of nodes has been selected for both dynamics.

$\bullet$ The state of the $2$-FS and $5$-FS is defined by: $f_2^{j,k}(t, u_j, w_k)$, $f_5^{j,k}(t, u_j, w_k)$, with $1<j<m$. $f_2^{j,k}$ and $f_5^{j,k}$ can be further characterized into $f_2^{j = 2,k}$, $f_5^{j = 2,k}$ and $f_2^{j\geq3,k}$, $f_5^{j\geq3,k}$ denoting, respectively, asymptomatic and symptomatic infection, while symptomatic infected, corresponding to higher values of $j$, might require hospitalization.

$\bullet$ The state of the $1, 3, 4$-FSs is defined by: $f_1^{1,k}(t, u_1, w_k)$, $f_3 = f_3(t)$ for individuals who succeed to reach back to the state $j = 1$; $f_4 = f_4(t)$ for infected individuals who, by reaching the state $j = m$, do not successfully recover.

Remark 2.1. If birth events are not considered, the normalization, with respect to $N_0$, preserves, for all $t \geq 0$ conservation of the number of individuals:

Remark 2.2. It is useful introducing the density of each FS defined as the number of individuals denoted, for each $i$-FS, by $n_i = n_i(t)$:

while $n_{3,4}(t) = f_{3,4}(t)$.

The dynamics of transition across FSs is shown in the flow chart of Fig. 1, where each block identifies the specific FS.

2.2.   On the derivation of models

Let us now transfer the aforementioned qualitative assumptions into a differential system according to detailed assumptions on the interactions involving a-particles of all FSs. Interactions modify the internal activity variables and, in addition, promote transition across FSs.

The kinetic theory of active particle with discrete states provides mathematical tools to model the dynamics of the dependent variables $f_i^{j,k} = f^{j,k}_i(t)$, where the latter refers to the probability to find an a-particle of the $i$-FS in the state $j,k$. The dynamics of interactions can be described by the following terms:

$\bullet$ $\eta_{pr}^{qs}(y,z)$ which denotes the interaction rate between the $y$-FS and $z$-FS, with microstates $p,r$ and $q,s$, respectively. Interactions promote transition from the $y$-FS to the $i$-FS.

$\bullet$ $\mathcal{A}_{pry}^{qsz}(y \to i, \,pr \to jk)$ which denotes the transition probability into the $i$-FS with microstate $j,k$ from the $y$-FS with microstate $p,r$ due to interactions with the $z$-FS with microstate $q,s$.

The formal structure that provides the conceptual framework to derive models that can be obtained by the superposition of kinetic type model, where the interaction rate and the transition probability determine the dynamics of the dependent variable, and a source term modeling transition across FSs and external actions, for instance vaccines:

where $\mathcal{G}^{j,k}_i[ \mathbf{f}]$ and $\mathcal{L}^{j,k}_i[ \mathbf{f}]$ denote the gain and loss terms of particles which, at time $t$, gain the state $j,k$ and which lose such a state, respectively, due to interactions across FSs, while $\mathcal{S}^{j,k}_i[ \mathbf{f}]$ denotes a source term modeling the inlet/outlet related to the $i$-FS due to the dynamics across FSs.

The kinetic theory of active particles provides the following expressions gain and loss terms:

and

while the modeling of the source term $\mathcal{S}$ depends on the transitions across functional subsystems. A phenomenological interpretation of biological reality can lead to a detailed modeling of the various interaction terms to be substituted into Eq. (3) thus obtaining specific models. In more details, we derive a model based on the following assumptions:

1. Interactions between $1$-FS and $2$-FS or $5$-FS and infection dynamics: Active $1$-FS particles interacting with a-particles from $2$-FS and $5$-FS might become, in probability, infected. The rate of infection depends on the level of progression $u_j$ of the infected individuals by a probability of infection which grows with $u_j$.

2. Interaction rate and social distance: The interaction rate depends on the local density somehow related to a social distance depending on time. Interactions do not modify the immune defence ability, while particles which move from $1$-FS to $2$-FS take the value $u_2$, corresponding to low level of infection.

3. In-host immune competition: Virus particles progress (proliferate) thanks to foraging of the surrounding tissues of the lung, while the immune defence counteracts the progression by inducing a regression.

4. Onset of variants: The derivation of models can be developed by adding a new population, namely a new FS, and which can generate infected individuals with higher infectivity ability.

5. Transition across FSs: A-particles from $2$-FS and $5$-FS move to $3$-FS if the immune defence succeeds to obtain a regression down to $u_1$, while a-particles from $2$-FS and $5$-FS move to $4$-FS if the immune defence does not succeed to obtain a regression. Individuals who succeed in recovering are supposed not to be susceptible to a new infection.

The description in Items 1–5 can be transferred into a mathematical framework by a detailed modeling of interactions which inserted them into the general structure (3) leading to models of virus dynamics. In more details, we consider the following models of interactions:

$\bullet$ Infection dynamics: Healthy individuals of the $1$-FS, with state $u_1$, interact with an infected individual from $2$-FS with state $u_j, \, j >1$, and becomes infected with a probability which depends on a parameter $\alpha = \alpha (t)$ and on the level of infection of the individual from $2$-FS.

$\bullet$ Interaction rate: The interaction rate is modeled by $\eta = \eta_0(\rho) \, \alpha (t)$, where $\eta_0$ depends on the local density of individuals $\rho$ (it grows with increasing value of the density thus reaching the maximal value in the case of maximal packing density).

$\bullet$ Loss of healthy people due to contagion: The dynamics refers to $f_1^{1,k}(t)$, interactions generate the following loss term:

for $k = 1, \ldots, m$, while we have put $\eta_0 = 1$ that is inserted in the scaling of the time variable. Analogous calculations, which are not repeated here, refer to the interactions between $1$-FS and $5$-FS.

$\bullet$ In-host dynamics of infected individuals: Each infected individual is a carrier of a struggle between virus particles and immune system. The virus takes advantage from the foraging by surrounding tissues and increases its micro-state from each $j$-level to the higher $(j+1)$-level depending on a parameter $\beta$ and on the $j$-th level. The immune system acts to decrease the $j$-level to the lower $(j-1)$-level depending on a parameter $\gamma_k$ and on the $k$-th level. Individuals, whose virus progression levels reach the values $u_1$ and $u_m$ move, respectively, to $3$-FS and $4$-FS.

The dynamics refers to $f_2^{j,k}(t)$ and it is governed by both gain and loss terms:

Two positive defined parameters $\beta \in [0,1]$ and $\gamma_k \in [0,1]$, where $\gamma_k$ depends on each $k$ level, are introduced to model, respectively, the proliferative ability of the virus and the defence ability of the immune system accounting also for the in-host interaction rate. Detailed calculations yield:

and

for $j = 2,\ldots,m-1$, $k = 1, \ldots, m$, and $\delta$ denotes the Dirac delta function.

$\bullet$ Dynamics of new variants: The specific features of the new FS are labeled by the subscript $5$ as follows:

where $\lambda$ to be related to the distance, measured by an appropriate metric in a vector space, at the molecular scale between the original virus and its variant. The dynamics can be described by the same approach of the primary virus, i.e. $\partial_t f_5^{j,k}(t) = \mathcal{G}_{5}^{j,k}[ \mathbf{f}](t) - \mathcal{L}_{5}^{j,k}[ \mathbf{f}](t)$, where the gain and loss terms are computed as (8) and (9) simply replacing $f_2$ by $f_5$ and $\beta$ by $\beta(1 + \lambda)$.

$\bullet$ Trend to recovering: The dynamics consider the inflow of individuals from $1$-FS and $5$-FS into $3$-FS:

$\bullet$ Need of hospitalization: The need of hospitalization is given by

for high values of the level of progression of the virus, i.e. increasing values of $j$ denote increasing values of the pathology.

$\bullet$ Trend to death: The dynamics are caused by the inflow from $2$-FS and $5$-FS into $4$-FS:

Collecting all equations into a differential system yields:

where $j = 2,\ldots,m-1$ and $k = 1, \ldots, m$.

The mathematical model (12) accounts for mutations and selection into variants, which act an important role in the spread and control of the pandemics. We are interested in understanding how these two events modify the dynamics in absence of them which is as follows:

Remark 2.3. Equations $(12)$ and $(13)$ describe the dynamics of heterogeneous FSs. The model can also predict the dynamics of individuals who need hospitalization with respect to those whom home care might be sufficient to reach full recovery. The superscript k considers different levels of the immune strength, that can be related - for instance - to the range of age of individuals in the population. An approximation of reality can be obtained by supposing that the immune strength is the same for all individuals which corresponds to $\gamma \cong \gamma_k$.

Remark 2.4. In absence of virus infection or total lack of awareness of the risk of infection, the locking parameter, which may also be called social distancing, is taken $\alpha = 1$. Once the awareness of the risk of contagion appears then $\alpha$ is reduced to $\alpha_\ell < \alpha$ by a "locking" action applied both by crisis managers and individual awareness. After locking, the social distance might be decreased as modeled by a parameter $\alpha_d > \alpha_\ell$ with $\alpha_d <1$.

3.   Simulations in absence of variants

This section presents some simulations selected with the aim of understanding how the dynamics of the epidemics develops after a locking action. We consider the dynamics in absence of mutations corresponding to model (13). In more details, we consider, for different values of the locking parameter $\alpha_\ell$, firstly the study of the shape of a first wave generated by an initial growth of the number of infected individuals which is followed by decay related to locking. Subsequently, we investigate the onset and shape of a second wave generated by a locking relaxation after the decay of the first wave. We study the shape depending on the locking parameter $\alpha_d$. Simulations indicate that there exists a critical value $\alpha_c$ such that if $\alpha_d > \alpha_c$ the height of the second wave is higher than that of the first wave, while the opposite appears when $\alpha_d < \alpha_c$.

All simulations are developed for $m = 11$, while the time variable is normalized with respect to the time $T^*$ corresponding to full decay of the first wave corresponding to the first simulation. We consider a primary virus which at initial time $t = 0$ is identified with an initial condition $\varepsilon$ small with respect to one, i.e. a small number of infected among the whole populations. Simulations refer to different values of the ratio $\kappa = \frac{\beta}{\gamma}$ and of the locking parameters $\alpha_\ell$ and $\alpha_d$, respectively. The preliminary study in this section, completed with that of the next section, can contribute to a detailed parameters sensitivity analysis which, at present, is still in progress.

$\bullet$ Let us consider a locking dynamics applied at $t = 0$ following the onset of a virus infection described by model (13) with initial conditions:

Figure 2 studies the influence of $\alpha_\ell$ and shows time dynamics of the infected population for $\alpha_\ell = 0.2, 0.3$, and $\kappa = 0.16$, while Figure 3 studies the influence of $\beta$ corresponding to $\alpha_\ell = 0.3$, and $\kappa = 0.1, 0.2$. As it is expected, increasing the locking actions, i.e. decreasing $\alpha_\ell$, decreases the level of infection, while increasing the proliferative ability of the virus, i.e. increasing $\kappa$, increases the level of infection. The slope of the increasing dynamic is sharp, while the decay is smooth. i.e. there is a loss of symmetry.

$\bullet$ Consider now the study of the onset of a second wave that is generated by a relaxation of the locking parameter. We refer to model(13) with initial conditions (14), with $\varepsilon = 0.001$, while the dynamical response is studied for $\kappa = 0.1$, de-locking time $T_d = 1$, and with $\alpha$ depending on time as follows:

while different values of $\alpha_d$ are selected.

Sample simulation are reported in Figs. 4, 5 which report the dynamic of $n_2$ versus time and show how a second wave appears when the locking parameter $\alpha_\ell$ is substituted by $\alpha_d > \alpha_\ell$. In more details, simulation shows how the second wave can rapidly reach levels worse than the first one, by a peak whose top increases with $\alpha_d$, see Fig. 4. On the contrary, Fig. 5 shows how the second wave can be lower than the first one if the locking is more restrictive, namely limited to $\alpha_d = 0.20, 0.25. 0,30$.

These simulations show an emergent behavior that has been effectively observed in the real situation, namely that although a locking is applied and the virus is brought to very low levels, apparently it has disappeared, but it remains at a latent state ready to start again once the social confinement is relaxed. This emergent behavior is confirmed by the onset of subsequent waves which always appear when a de-locking action is applied. Therefore, an important indication to crisis managers is that the locking parameter can be brought to a value $\alpha_d > \alpha_\ell$, but an excessive increase of $\alpha_d$ brings to a second wave much worse than the first one. An additional feature shown by these simulations in the case of high values of $\alpha_d$ the second wave rapidly appears, while decreasing values of $\alpha_d$ delay the onset of the second wave. Additional simulations, which are not reported here, confirm what was already observed in [7]. Namely if de-locking is applied during the descending phase of the first wave before having reached minimum value of the infected people, then a new wave rapidly appears with a peak that increases with increasing $\alpha_d$. This behavior is studied in the next section for a dynamics which includes the onset of variants of the primary virus.

4.   Simulation with mutations dynamics

This section presents a selection of simulations which aim at understanding how the dynamics of the epidemics develops in the presence of mutations. These generate the onset of variants, whose dynamic is described by model (12), that modifies the dynamics of the primary virus described by model (13). Simulations are referred to the locking strategy with the objective of showing, by quantitative results, how the more aggressive variants end up with replacing the primary virus. As in Section 3, simulations are developed for $m = 11$ referring to the locking parameters $\alpha_\ell$ and $\alpha_d$, and to the parameter $\kappa$, while the initial conditions include the initial state of the variant $\varepsilon_v$.

Empirical evidence shows that mutations can generate aggressive variants that show an ability to infect greater than that of the primary virus. This general trend can be explained by the model proposed in our paper as it is based on the assumption that a virus with high proliferative ability produces a high viral charge which increases the infectivity ability. Simulations show how far the variant increases the number of infected people, as well as the death of those who do not succeed in recovering.

$\bullet$ The dynamics of an early de-locking after the first wave is shown in Fig. 6, where the dynamics of $n_2$ and $n_5$ are reported versus time for locking time $T_d = 0.5$. Simulations depict how the density of the primary virus is still at a high value when the de-locking is applied hence it undergoes a high density second wave higher than that of the variant. However, the variant becomes, when time increases, equivalent, or even higher, than the primary. Figures 7 and 8 study the same dynamics of $n_2$ and $n_5$ corresponding to different values of $T_d$. In details, Fig. 7 shows how for $T_d \cong 1$, the predominance of the variant is immediate, while Fig. 8 shows an intermediate situation corresponding to the de-locking time $T_d \cong 0.75$, where the trend of the variant to dominate and replace the primary virus is already enhanced with respect to Fig. 6.

$\bullet$ Figures 9 and 10 show $n_4 = n_4(t)$ which defines the increase versus of individuals who are not able to contrast the progression of the virus and die. Figure 9 shows, in the case of absence of mutations, how $n_4$ grows during the first and second wave. Simulations are developed at fixed locking time $T_d = 1$, $\varepsilon = 0.001$ and $\kappa = 0.1$ and indicate how the number of dead people increases with increasing de-locking action for different levels of $\alpha_d = 0.4, 0.5, 0.6$. Figure 10 is somehow analogous, but it refers to the case of presence of mutations. Simulations can be referred to those of Figure 8 as the same values $T_d = 1$, $\varepsilon = 0.01$ and $\kappa = 0.1$ are used, while $\alpha_d = 0.5$. The model shows that the presence of a variant anticipates the onset of the second wave and increases the number of deaths. This number increases with the aggressiveness of the variant which is modeled by the parameter $\lambda$.

Simulations presented in this section do not claim that an exhaustive description has been given. Indeed, further simulations accounting for a detailed parameters sensitivity analysis should be developed. However, all simulations clearly show that the onset of a new variant, more aggressive than the primary virus, generates a progressive prevalence of the variant over the firstly appeared virus. This overcoming appears by different types of dynamics, but ultimately confirming the aforementioned trend. This is an important piece of information, as it contributes to the amount of information necessary to design strategies to fight the virus. Further reasonings follow in the critical analysis proposed in the next section.

5.   Critical analysis

The mathematical model proposed in this paper considers the multiscale dynamics of the virus from the infection of individuals to the within-host competition which may lead to recovery, but also to hospitalization, or even death. In more details, our model includes virus mutations that generate, by selection, new more aggressive variants and a multiscale modeling of the in-host vaccination dynamics related to strategies to mitigate the damage of pandemics [37].

Without repeating concepts already stated in the first sections, we stress that a multiscale approach is necessary to consider the aforementioned mutations and selections. Indeed, more aggressive variants of the virus exhibit higher proliferative ability which increases the infectivity properties of the virus. These dynamics are related to the fact that increasing the number of virus particles in the lung, namely the so-called viral load, makes the individual more infective when in contact with other individuals. Unfortunately, these dynamics end up with increasing the number of infected humans who need hospitalization or even die.

Comparing the present simulations to those developed for the simpler model [7], we have used a higher number of collocation nodes to account in a more precise way for the in-host dynamics and, specifically, the progression of the pathological state. Computing the number of individuals that need hospitalization is an important information to be delivered to crisis managers as in the pandemic time hospitals have been obliged to face crisis situations related to lack of space for patients affected by pathologies different from that of Covid19. Hence a higher number of internal states can refine the indications by referring the number of people who need hospitalization to their level of pathology.

The computational study has also allowed the study of emerging behaviors such as the onset of new waves and their shape depending on the parameters of the model. In addition, hospitalization dynamics can be related to different levels of progression of the pathology.

Therefore, the model appears to be able to provide a panorama of emerging behaviors which may contribute to the activity of crisis managers. The main focus of our simulations has been referring to the role of variants on the overall response of the system. In particular, the following specific features of the model, as well as of the plan of simulations, can be selected as appropriate to contribute to support crisis management:

$\bullet$ Heterogeneity of the population corresponding to the label $k$ with $k = 1, \ldots, m$ which can be related to different abilities of the immune defence. This feature can be referred to different levels of the age of the populations corresponding to levels of the ability to develop the immune defence against the virus.

$\bullet$ The study of locking and de-locking strategy can contribute to plan the applicable strategies, i.e. the levels $\alpha_\ell$ and $\alpha_d$, as well as the behavior versus time, by which locking and de-locking are applied in order to reduce the spread of the epidemics. In addition the study can be further specialized on the number of individuals who need hospitalization and on their trend to death.

$\bullet$ Understanding how the onset of variants adds to the original virus. Simulations have shown how variants can progressively take the place of the original virus and produce an increase of hospitalization flux and death.

However, we do not naively claim that the sample simulations we have presented in our paper do cover an exhaustive parameter sensitivity analysis. Simply, we have shown that the model can provide a rich description of different scenarios that might hopefully contribute to strategies that make lighter crisis situations related to the pandemics. A program to develop a systematic sensitivity analysis is in progress in order to design a database, where the broad variety of scenarios of dynamics of the complex system can be tested.

Hopefully, new studies might take advantage of the indications given above and contribute to further developing the model. The following specific modeling perspectives are selected among various possible ones:

(i) Heterogeneity referred not only to the ability of the immune defence, that can be somehow referred to age, but also to the role of individuals in society. This means accounting for both levels of contagion risk during work and of the economical impact on society.

(ii) Contagion in crowds starting from the pioneering contributions [12,34,35] towards development of models of crowd dynamics which include social awareness [4,9,43], where one of the key problem consists of modeling the contagion parameter [25]. Contagion in crowds described by a macroscopic equation derived from the underlying description delivered by kinetic theory [6].

(iii) Modeling vaccination to improve the ability of immune defence and medical actions to weaken the proliferative and invasive ability of the virus. The modeling approach should start from a knowledge of vaccine actions, for instance from [16,29,37,45] to the mathematical description of the dynamics within vaccinated individuals and viruses.

(iv) Modeling in-host space dynamics by reaction-diffusion equations [10], based on micro-macro derivation [15].

(v) Calibration of models based on empirical data referred to well defined regional areas. Indeed, empirical data have reached a sufficient level of precision and, thanks to this increase of reliability, an effective toning of SEIR appeared related to control problems [39].

Further developments of the model proposed in our paper can contribute to tackle the research perspectives that have been indicated above. These, as mentioned, represent a selection based on our personal ideas and research past experience. A common aspect of these perspectives is that a key dynamics, which is present in all of them, is the dynamics of the immune competition which has been included in our approach at a very primitive level. This topic definitely deserves further developments and improvements.

Acknowledgments

Nicola Bellomo acknowledges the support of the University of Granada, Project Modeling in Nature MNat from micro to macro. This paper has been partially supported by the Consejeria de Economia, Conocimiento, Empresas y Universidad and European Regional Development Fund (ERDF), ref. SOMM17/6109/UGR. https://www.modelingnature.org.