key: cord-0634316-r3l7ag77
authors: Auricchio, Ferdinando; Colli, Pierluigi; Gilardi, Gianni; Reali, Alessandro; Rocca, Elisabetta
title: Well-posedness for a diffusion-reaction compartmental model simulating the spread of COVID-19
date: 2022-03-21
journal: nan
DOI: nan
sha: 78f70ce99d18c51ecc469255b3f7610ff8a5e6e6
doc_id: 634316
cord_uid: r3l7ag77

This paper is concerned with the well-posedness of a diffusion-reaction system for a Susceptible-Exposed-Infected-Recovered (SEIR) mathematical model. This model is written in terms of four nonlinear partial differential equations with nonlinear diffusions, depending on the total amount of the SEIR populations. The model aims at describing the spatio-temporal spread of the COVID-19 pandemic and is a variation of the one recently introduced, discussed and tested in [A. Viguerie et al, Diffusion-reaction compartmental models formulated in a continuum mechanics framework: application to COVID-19, mathematical analysis, and numerical study, Comput. Mech. 66 (2020) 1131-1152]. Here, we deal with the mathematical analysis of the resulting Cauchy--Neumann problem: the existence of solutions is proved in a rather general setting and a suitable time discretization procedure is employed. It is worth mentioning that the uniform boundedness of the discrete solution is shown by carefully exploiting the structure of the system. Uniform estimates and passage to the limit with respect to the time step allow to complete the existence proof. Then, two uniqueness theorems are offered, one in the case of a constant diffusion coefficient and the other for more regular data, in combination with a regularity result for the solutions.

The COVID-19 outbreak that deeply affected the world starting from the first months of 2020 led to a new, strong interest by researchers towards the development of mathematical models of infectious diseases, allowing the assessment of different scenarios and aiming at assisting the complex political decision-making process during the pandemic. As a consequence, many papers were recently published, proposing interesting modeling ideas (see, e.g., [1, 9, 12, 13, 17, 18, 20, 23, 25] ), often based on compartmental models, where the considered population is divided into "compartments" based on their qualitative characteristics (like, e.g., "susceptible", "infected", "recovered"), with different assumptions about the nature and rate of transfer across compartments. Despite this kind of models do not readily offer the possibility of a multiscale vision (as, e.g., proposed in [3] ), that would be a preferred feature given the nature of the phenomena to be simulated, they have the advantage of allowing a relatively easy introduction of diffusion terms (in this context, compelling ideas on diffusion models can be found, among others, in [4, 5, 24] and references therein). For a recent overview of mathematical models for virus pandemic, interested readers are referred to the articles included in [6] . Here, a compartmental model able to describe the spatial heterogeneity underlying the spread of an epidemic in a realistic city network is proposed in [8] . Another recent example where a classical compartmental model is enhanced to take into account infected individuals traveling on lines of fast diffusion is instead given by [7] . In fact, compartmental models are typically governed by a system of ordinary differential equations (ODEs) in time, which might possibly be enriched with coupling terms, additional equations, or other approaches to take into account also the spatial variation of the studied problem. An interesting alternative is to resort to compartmental models based on partial differential equations (PDEs), to accurately account for spatial variations at the continuum level.

Following this latter idea, in the present contribution, we consider the PDE-based model introduced in [21, 22] (a variation of it with delay terms is discussed in [16] ). In particular, the following system has been considered ∂ t s = αn − (1 − A 0 /n)β i si − (1 − A 0 /n)β e se − µs + div(nν s ∇s) (1.1) ∂ t e = (1 − A 0 /n)β i si + (1 − A 0 /n)β e se − σe − φ e e − µe + div(nν e ∇e) (1.2) ∂ t i = σe − φ d i − φ r i − µi + div(nν i ∇i) (1.3) ∂ t r = φ r i + φ e e − µr + div(nν e ∇e) (1.4) ∂ t d = φ d i .

(1.5)

In the above equations, the symbols s, e, i, r, and d denote the susceptible population, the exposed population, the infected population, the recovered population, and the deceased population, respectively, and n := s + e + i + r is the living population. Note that d refers only to deaths due to Due to the names of the compartments used, this model may be called a susceptible-exposed-infected-recovered-deceased (SEIRD) model. Moreover, equations (1.1)-(1.4) are complemented by initial conditions and no-flux boundary conditions, while just an initial condition is associated to (1.5) .

In [21] the authors presented a SEIRD mathematical model based on partial differential equations coupled with a heterogeneous diffusion model. The model was used to describe the spatio-temporal spread of the COVID-19 pandemic and to capture dynamics also based on human habits and geographical features. To test the model, the outputs generated by a finite-element solver were compared with measured data over the Italian region of Lombardy and the results showed a strong qualitative agreement between the simulated forecast of the spatio-temporal COVID-19 spread in Lombardy and epidemiological data collected at the municipality level.

In [22] the authors proposed a formulation of compartmental models based on partial differential equations through familiar continuum mechanics concepts, interpreting such models in terms of fundamental equations of balance and compatibility, joined by a constitutive relation. Such an interpretation was useful to aid understanding and interdisciplinary collaboration. The authors formally derived the model sensitivity to diffusion, described its growth and decay, and established its stability in the L 1 norm. Attention was paid to an ODE version of the model, using it to derive a basic reproduction number R 0 as well as analyzing its spectrum. Additionally, a series of numerical simulations showed the role that numerical methods, diffusion, and specific ingredients played in the behavior of the system. Implicit models were found to be effective in describing the temporal dynamics of the system, and second-order in-time methods in particular. The model proposed in [22] and used for simulating the COVID-19 spread in Lombardy, was further exploited in [14, 15] . In particular, finite element simulations were carried out in the context of adaptive mesh refinement and coarsening, forecasting the COVID-19 spread also in the U.S. state of Georgia and in the Brazilian state of Rio de Janeiro. Good agreements with real-world epidemiological data were attained in both time and space. Coming back to system (1.1)-(1.5), we note that d can be determined from (1.5) and the associated initial condition whenever the system of the first four equation in solved. Thus, since we are just interested in well-posedness, we do not consider the last equation and study system (1.1)-(1.4). Namely, we consider a slightly modified system. Indeed, we replace the terms 1 − A 0 /n and φ d i appearing in (1.1)-(1.2) and (1.3) by a non-singular one of the form A(n) and by the product φ d in, respectively. On the contrary, we can treat the degeneracy in the diffusion terms that occur as n approaches zero, by proving that n is bounded away from zero. This is done even in the more general situation of a coefficient of the form κ(n), where κ is continuous and strictly positive on (0, +∞). Finally, we assume that the constants ν s , ν e , ν i and ν r are the same. For this modified model, we prove an existence result under mild assumption on the initial data. Furthermore, we give two uniqueness results: the main assumption for the first one is that κ is a positive constant, while the second one requires that the nonlinearities and the initial data are smooth. Under the latter assumptions, we prove a regularity result as well.

We hope that our results could capture the interest for the resulting systems and give rise to new approximations and simulations, in order to deal with and test our variations in the framework of some geographical areas. Moreover, we claim that our arguments for well-posedness provide a serious mathematical validation to this class of SEIR models.

The paper is organized as follows. In the next section, we list our assumptions and notations and state our results. The proof of Theorem 2.2 regarding the existence of a solution is given in Section 4 and is prepared in Section 3 where an approximating problem obtained by time discretization is studied. Finally, Section 5 is devoted to the uniqueness of the solution.

In this section, we state precise assumptions and notations and present our results. First of all, the subset Ω ⊂ R 3 (lower-dimensional cases can be treated in the same way) is assumed to be bounded, connected and smooth. The symbol ∂ ν denotes for the normal derivative on Γ := ∂Ω. Moreover, we set for brevity Q t := Ω × (0, t) and Q := Ω × (0, T ).

(2.1)

If X is a Banach space, · X denotes both its norm and the norm of X d . The only exception from this convention on the norms is given by the spaces L p (1 ≤ p ≤ ∞) constructed on (0, T ), Ω and Q, whose norms are often denoted by · p , and by the space H defined below and its powers, whose norms are simply denoted by · . We put Moreover, V * is the dual space of V and · , · is the dual pairing between V * and V .

In the sequel, we work in the framework of the Hilbert triplet (V, H, V * ) obtained by identifying H with a subspace of V * in the usual way. Thus, by also using the symbol ( · , · ) for the standard inner product of H (this symbol will be used for any power of H as well later on), we have g, v = (g, v) for every g ∈ H and v ∈ V . Now, we list our assumptions on the structure of our system. We just consider the functions A and κ mentioned in the introduction and the constants α and µ, since all the other constants appearing in (1.1)-(1.4) will be normalized in the following. However, we keep them for a while, for the reader's convenience. We assume that A, κ : (0, +∞) → R are continuous, with A nonnegative and κ strictly positive (2.3) α and µ are positive constants (2.4) and we observe that a term like κ(n) = n is allowed in the equations. In principle, the system we are interested in is the following ∂ t s + A(n)β i si + A(n)β e se + µs − div(κ(n)∇s) = αn (2.5) ∂ t e − A(n)β i si − A(n)β e se + σe + φ e e + µe − div(κ(n)∇e) = 0 (2.6)

where each equation is complemented by no-flux boundary conditions and an initial condition. However, it is more convenient to consider an equivalent system in the unknowns n, s, i and h, where n and h are related to the other unknowns by n := s + e + i + r and h := s + e .

(2.9)

The new system is obtained from the previous one by summing up all the equations, keeping (2.5) and (2.7) and adding (2.6) to (2.5) . Hence, it is given by

(2.13)

However, as said above, among all of the constants appearing in these equations, just α and µ will play a significant role in the mathematical treatment. Hence, we replace the other constants and even some sums of them by 1, without loss of generality. At this point, we are ready to give our notion of solution. For the initial data we require that

and we call solution a quadruplet (n, s, i, h) enjoying the requirements

and satisfying, a.e. in (0, T ) and for every v ∈ V , the variational equations

and the initial conditions So we do not mind about it in the following. We also observe that the above variational equations are equivalent to their integrated versions with time dependent test functions. For instance, (2.17) is equivalent to

Here is our first result. As for uniqueness, we present two results. The first one regards a particular case and its proof, given in Section 5.1, is elementary. Our last result regards uniqueness in the case of a non-constant κ. Its proof, given in Section 5.2, is much more involved and is strictly related to a high regularity of the solution. For this reason, we have to assume that both the nonlinearities and the initial data are much smoother.

A is locally Lipschitz continuous and κ is a C 1 function

Then, the solution to problem (2.17)-(2.21) satisfying (2.15)-(2.16) is unique and enjoys the following regularity properties n, s, i, h ∈ W 1,p (0, T ; L p (Ω)) ∩ L p (0, T ; W 2,p (Ω)) for every p ∈ [1, +∞).

(2.25)

Throughout the paper, we make use of the Hölder inequality and the Sobolev inequality related to the continuous embedding V ⊂ L p (Ω) with p ∈ [1, 6] (since Ω is three-dimensional bounded and smooth). We also account for the elementary identity and inequalities

and quote (2.27) as the Young inequality. Furthermore, we take advantage of the summation by parts formula

which is valid for arbitrary real numbers a 1 , . . . , a m and b 0 , . . . , b m .

In this section, we prepare the proof of Theorem 2.2 by introducing and solving an approximating problem obtained by time discretization. However, the structural functions A and κ have to satisfy different assumptions and the initial data have to be smoother.

In the next section, by starting from the original structure and the original initial data, we consider the discrete problem with structural functions and approximating initial data constructed in order to satisfy the assumptions listed below. Two constants κ * and κ * and two real functions A and κ defined in the whole of R are given such that

Notation 3.1. However, we prefer to use the lighter symbols A and κ instead of the heavy A and κ. Indeed, no confusion can arise since the original functions A and κ introduced in (2.3) will never appear within the section.

For a fixed positive integer N, we set τ := T /N. Then, the time-discretized problem we are going to study is the following: given

we look for four (N + 1)-tuples

satisfying, for k = 0, . . . , N − 1, the variational equations

all for every v ∈ V , as well as the sign and initial conditions

Theorem 3.2. Under assumptions (3.1) and (2.4) on the structure, suppose that

Then, if the initial data satisfy (3.2), problem (3.4)-(3.9) has a unique solution.

We prepare an easy lemma. Then, the problem of finding u ∈ V satisfying the variational problem

has a unique solution, and this solution satisfies

in Ω, then u ≥ λ a.e. in Ω.

Proof. The existence of a unique solution is clear since the left-hand side of (3.12) is an inner product in V that is equivalent to the usual one. Moreover, the first inequality of (3.13) is given by the weak maximum principle and the same can be said for the last sentence of the statement, since the equation obtained by replacing f by f −λ b is satisfied by u − λ . Now, we prove the second inequality of (3.13). We set w :

we have w ≤ 0 by the weak maximum principle, i.e., the desired inequality.

Proof of Theorem 3.2. It suffices to prove the following: for k = 0, . . . , N − 1, if (n k , s k , i k , h k ) belongs to (V ∩ L ∞ (Ω)) 4 and satisfies the inequalities in (3.8), then system (3.4)-(3.7) has a unique solution (n k+1 , s k+1 , i k+1 , h k+1 ) belonging to (V ∩ L ∞ (Ω)) 4 and satisfying the analogous inequalities, i.e., n k+1 ≥ 0, h k+1 ≥ s k+1 ≥ 0 and i k+1 ≥ 0. We recall (3.2) for the case k = 0, fix k and (n k , s k , i k , h k ) as said and show that we can find a unique solution to (3.4)-(3.7) with the proper sign conditions. This is done in two steps.

Solution to the first equation. We introduce the function K : R → R by setting

and we would like to assume u := K(n k+1 ) as the new unknown for (3.4). Due to (3.1) for κ, the function K is one-to-one, onto and Lipschitz continuous, and its inverse is Lipschitz continuous too. Thus, equation (3.4) can be rewritten in term of u. Namely, noting that ∇u = ∇(K(n k+1 )) = κ(n k+1 )∇n k+1 , we can write it as the variational formulation of the homogeneous Neumann problem for the equation

Notice that every variational solution automatically belongs to W and satisfies both (3.15) and the homogeneous Neumann condition due to elliptic regularity. We claim that the new problem has a unique solution. Indeed, it is equivalent to the minimization of the functional J : V → R given by

and we show that this problem has a unique solution. Notice that J actually is welldefined, since K −1 is Lipschitz continuous, λ ∈ L ∞ (Ω) and n k ∈ H. Moreover, J is strictly convex and coercive. To see this, it suffices to recall that inf λ > 0 (since we are assuming (3.10) and i k is nonnegative) and observe that

Hence, the above problem has a unique solution u and we conclude that n k+1 = K(u) is the unique solution to (3.4) . Since u ∈ W , we deduce that u is bounded so that that n k+1 is bounded too. To see that it is nonnegative, it suffices to rearrange (3.4) and apply Lemma 3.3.

Solution to the other equations and conclusion. Since n k+1 is known, we can solve (3.5) for s k+1 by applying the first part of Lemma 3.3. We can do the same to find h k+1 and i k+1 from (3.7) and (3.6) in this order. The lemma also ensures that s k+1 , h k+1 and i k+1 are bounded and nonnegative, provided we can prove that h k+1 ≥ s k+1 . To this end, we set e k := h k − s k and e k+1 = h k+1 − s k+1 , take the difference between (3.7) and (3.5) and write it in the form 

This section is devoted to the proof of Theorem 2.2. Our argument relies on a priori estimates on the solution to a suitably specified discrete problem and on the convergence of proper interpolating functions. So, we introduce some notations concerning interpolation at once. 

For the reader's convenience, we summarize the relations between the finite set of values and the interpolating functions in the following proposition, whose proof follows from straightforward computation:

Moreover, it holds that

and similar identities for the difference z τ − z τ . Finally, we have that

if z ∈ H 1 (0, T ; Z) and z k = z(kτ ) for k = 0, . . . , N.

(4.11)

We are now ready to properly specify the discrete problem. Namely, starting from A and κ as in (2.3), we introduce new functions A and κ obtained by a truncation operator T to be used in the discrete problem. We set n * := e 2T (α−µ) + n 0 ∞ , s * := s 0 ∞ + T αn * h * := h 0 ∞ + T (αn * + s * ) i * := i 0 ∞ + T (h * + s * ) and n * := e −T (i * +(µ−α) + ) inf n 0 (4.12)

and define A, κ : R → R by setting for y ∈ R A(y) = A( T (y)) and κ(y) = κ( T (y)) where T (y) := max{n * , min{y, n * }}. (4.13)

Next, we approximate the initial data n 0 , s 0 , i 0 and h 0 as in (2.14) by smoother functions n 0 τ , s 0 τ , i 0 τ and h 0 τ satisfying n 0 τ ∈ V ∩ L ∞ (Ω) and 0 ≤ n 0 τ ≤ n 0 ∞ a.e. in Ω (4.14) n 0 τ ≤ n 0 and τ n 0 τ 2 V ≤ n 0 2 (4.15) n 0 τ → n 0 strongly in H as τ ց 0 (4. 16) and the analogues for s 0 τ , i 0 τ and h 0 τ as well as h 0 τ ≥ s 0 τ and n 0 τ ≥ inf n 0 a.e. in Ω. (4.17) This can be done, e.g., by a singular perturbation argument. Indeed, if τ ∈ (0, 1) and u ∈ L ∞ (Ω) is nonnegative, the unique solution u τ ∈ V to the variational problem

V ≤ u 2 , and converges to u strongly in H as τ ց 0. Finally, we assume

The functions A and κ satisfy (3.1) with κ * := min{κ(y) : n * ≤ y ≤ n * } and κ * := max{κ(y) : n * ≤ y ≤ n * }. Proof. All the components n k are nonnegative. However, as A and κ are defined in the open interval (0, +∞), we have to reinforce this and show that each n k is strictly positive in order to give a meaning to A(n k ) and κ(n k ). We prove that n * ≤ n k ≤ n * a.e. in Ω for k = 0, . . . , N . Since A(y) = A(y) and κ(y) = κ(y) for every y ∈ [n * , n * ], this also yields that A(n k ) = A(n k ) and κ(n k ) = κ(n k ) for k = 0, . . . , N, i.e., the thesis of the statement.

The proof of (4.21) is done in several steps. The last of them needs L ∞ bounds for the other components of the solution, and these are proved as well.

First upper bound. We assume 0 ≤ k < N and check that we can apply Lemma 3.3 to n k+1 . Indeed, (3.4) with κ in place of κ implies that n k+1 solves (3.12) with

Moreover, we can take a 0 = κ * and b 0 = (1/τ ) − (α − µ) + since i k is nonnegative. Notice that b 0 > 0 by (4.18) . As n k+1 is nonnegative, the lemma implies that

As this holds for 0 ≤ k < N, n 0 = n 0 τ and n 0 τ ∞ ≤ n 0 ∞ by (4.14), we deduce that

for k = 0, . . . , N. Now, observe that ln(1 − y) ≥ −2y for every y ∈ (0, 1/2]. Indeed, the logarithm is a concave function and the inequality is satisfied at y = 1/2 and for y > 0 small enough by comparison of the derivatives at 0. Since τ (α − µ) + ≤ 1/2 by (4.18), it follows that

Therefore, the second inequality of (4.21) is proved. Now, we prove analogous upper bounds for the other components of the solution by applying Lemma 3.3 once more, always with a = κ(n k+1 ).

Further upper bounds. We start with s k . Equation 

By recalling that A is nonnegative and that i k ≥ 0 and h k ≥ s k a.e. in Ω, we can take b 0 = (1/τ ) + µ. Since s k+1 is nonnegative, the lemma yields

On account of (4.24) and (4.12), we deduce that s k ∞ ≤ s 0 ∞ + kτ αn * ≤ s * for k = 0, . . . , N. For h k+1 we can take b = 1 τ + 1 + µ and f = 1 τ h k + α n k+1 + s k+1 .

Then the lemma and the estimates for n k and s k already proved give

Hence

Finally, for i k+1 we can take b = 1 τ + n k+1 + 1 + µ and f = 1

and a quite similar calculation also using (4.26) yields

Lower bound. Finally, we prove the left inequality of (4.21). We fix k with 0 ≤ k < N and assume that n k ≥ λ k a.e. in Ω for some λ k > 0 and we prove that n k+1 ≥ λ k+1 a.e. in Ω where

.

As already observed, (3.12) is solved by n k+1 with a, b and f given by (4.22) . Let us estimate f − λ k+1 b from below. We have

Then, we can apply the last sentence of Lemma 3.3 and obtain that n k+1 ≥ λ k+1 a.e. in Ω.

Since n 0 = n 0 τ ≥ inf n 0 a.e. in Ω, we can take λ 0 = inf n 0 and have

a.e. in Ω, for k = 0, . . . , N.

(4.28)

Since ln(1 + y) ≤ y for every y > −1, we obtain

Hence, (4.28) yields that n k ≥ n * a.e. in Ω for every k. This concludes the proof.

At this point, we establish two estimates whose proof is made very easy by the above L ∞ -bounds. For the second one it is convenient to rewrite the estimates in terms of the interpolating functions. First a priori estimate. We rewrite (3.4)-(3.7) in the form n k+1 − n k τ , v + (κ(n k+1 )∇n k+1 , ∇v) + κ * (n k+1 , v) = (g n,k , v) (4.29)

all for every v ∈ V , where we have set

We observe that (4.21), (4.25)-(4.27) and the continuity of A imply that

Now, we test (4.29) by τ n k+1 and sum over k from 0 to m − 1, with an arbitrary positive integer m ≤ N. We obtain

(g n,k , n k+1 ).

Since κ ≥ κ * , by owing to (2.26) and the Schwarz and Young inequality, we deduce that

Then, it suffices to rearrange and account for the first condition in (4.15) and (4.33) to obtain max k=0,...,m 

At this point, we introduce the interpolating functions. Namely, we term n τ , s τ , i τ , h τ , n τ , s τ , i τ , h τ and n τ , s τ , i τ , h τ 

and the analogues for n τ , s τ , i τ and h τ . Next, the estimates (4.34)-(4.37) imply n τ L 2 (0,T ;V ) + s τ L 2 (0,T ;V ) + i τ L 2 (0,T ;V ) + h τ L 2 (0,T ;V ) ≤ c (4.40) n τ − n τ L 2 (0,T ;H) + s τ − s τ L 2 (0,T ;H)

Notice that (4.41) and (4.10) provide the same estimate for n τ − n τ . Moreover, the property (4.8) with Z = V , the second inequality in (4.15) and (4.39) also yield

In order to obtain an estimate for the time derivatives and let τ tend to zero, we write equations (3.4)-(3.7) in terms of the interpolating function. We have

all of them being satisfied for every v ∈ V a.e. in (0, T ). Moreover, requirements (3.8) and the initial conditions read Second a priori estimate. We take an arbitrary v ∈ L 2 (0, T ; V ), write (4.43) at the time t and test it by v(t). Then, we integrate over (0, T ) and rearrange. Thanks to the above estimates, we obtain

By analogously treating the other equations, we conclude that

Since all this holds under the only restriction on τ given by (4.18), we are ready to let τ tend to zero. as τ ց 0 (more precisely for a suitable sequence τ j ց 0). In particular, the values at 0 converge weakly in H to the left-hand sides of (2.21) and the initial conditions (2.21) themselves are satisfied on account of (4.16). Moreover, the differences n τ −n τ and n τ − n τ converge to zero strongly in L 2 (0, T ; H), and also a.e. in Q without loss of generality, so that n τ and n τ tend to n strongly in L 2 (0, T ; H). Since the same holds for the other variables, all the inequalities in (2.16) are satisfied too. Moreover, (4.49) yield n * ≤ n ≤ n * , s ≤ s * , i ≤ i * and h ≤ h * a.e. in Q.

By combining with the point-wise convergence, we infer that n τ , s τ , i τ and h τ converge to their limits strongly in L p (Q) for every p ∈ [1, +∞) and the same holds for n τ , s τ , i τ and h τ . We have some consequences. First, the products like i τ n τ in (4.43) converge to the right products (in this case the limit is in) in L p (Q) for p ∈ [1, +∞). Next, A(n τ ) converges to A(n) in the same topology, since it is bounded in L ∞ (Q) and converges to A(n) a.e. in Q, all this for A is continuous. By the same reason, κ(n τ ) converges to κ(n) strongly in L p (Q) for p < +∞. From this, we claim that κ(n τ )∇n τ → κ(n)∇n, κ(n τ )∇s τ → κ(n)∇s, κ(n τ )∇i τ → κ(n)∇i and κ(n τ )∇h τ → κ(n)∇h weakly in (L 2 (0, T ; H)) 3 .

In fact we prove the first one, only, since the others are analogous. For the product, we have weak convergence to κ(n)∇n in (L q (Q)) 3 for q ∈ [1, 2). On the other hand, κ(n τ ) is bounded in L ∞ (Q) and ∇n τ is bounded in (L 2 (0, T ; H)) 3 , so that the product is bounded in (L 2 (0, T ; H)) 3 and thus has a weak limit in this topology. Clearly, this weak limit has to be κ(n)∇n.

All this allows us to let τ tend to zero in (4.43)-(4.46). We consider the first one. We write its integrated version, namely

Then, it is clear that we can let τ tend to zero and obtain (2.22) . As the same argument works for the other equations, the proof is complete.

In this section we give the proofs of Theorems 2.3 and 2.4.

We assume that A is locally Lipschitz continuous and that κ is a positive constant. We pick any two solutions (n j , s j , i j , h j ), j = 1, 2, and prove that they are the same. First, we make some observations. By recalling that n j are bounded (as well as the other components) and that inf n j > 0, we fix an interval [n ⋆ , n ⋆ ] ⊂ (0, +∞) that contains all the values of n 1 and n 2 . Then, the restriction of A to [n ⋆ , n ⋆ ] is Lipschitz continuous. As for the estimates we are going to prove, we reinforce the convention given in Notation 4.4 for the constants by allowing the values of c to depend on the fixed solutions (through their L ∞ norms), in addition.

We set for brevity n := n 1 − n 2 , s := s 1 − s 2 , i := i 1 − i 2 and h := h 1 − h 2 . We write (2.17) for both solutions, take the difference and test it by n. Then, we integrate over (0, t) and adjust. By also owing to the Young inequality, we have

By proceeding with (2.18) in the same way, we obtain At this point, it suffices to add (5.1)-(5.3) to each other and apply the Gronwall lemma to conclude that n = s = i = h = 0.

It is understood that the assumptions listed in the statement are in force. Since the component n of every solution belongs to L ∞ (Ω) and is bounded away from zero (see (2.15 )-(2.16)), we can assume, without loss of generality, that κ * ≤ κ(y) ≤ κ * for some positive constants κ * and κ * and every y > 0 (5.4) and that A is Lipschitz continuous, whenever we fix one or two solutions. The bounds κ * and κ * and the Lipschitz constant of A depend on the solutions we fix every time, and the same holds for the Lipschitz constants of the functions we are going to introduce. We set and observe that the regularity and boundedness properties of u, but that of the time derivative, are the same as those of n.

Our project is the following: first we prove that every solution (n, s, i, h) to problem (2.17)-(2.21) satisfying (2.15)-(2.16) enjoys the regularity properties specified in (2.25); then, we prove uniqueness on account of the regularity already proved. The main difficulty is a regularity result for the component n of the solution. This needs some preliminary results. 

Proof. We fix a solution (n, s, i, h) and suitably extend the component n. We also introduce an auxiliary function. We define n and f on (−1, T ) by setting Notice that div(κ(n 0 )∇n 0 ) belongs to H by (2.24) (for this it would be sufficient a much weaker assumption). Also notice that n is a continuous H-valued function so that no jumps for its derivative can occur at t = 0. Therefore, n and f satisfy n ∈ H 1 (−1, T ; V * ) ∩ L 2 (−1, T ; V ) and f ∈ L 2 (−1, T ; H)

However, for simplicity, we write n in place of n in the sequel. Now, we fix τ > 0 small (namely, τ < min{1, T }). Later on, we let τ tend to zero. We set u(t) := K(n(t)) and U(t) := by noting that ∇u = κ(n)∇n and ∂ t U = u. Then, for t ∈ (−τ, T − τ ), we integrate (5.10) on (t, t + τ ) with respect to time and test the resulting equality by τ −2 (u(t + τ ) − u(t)).

We obtain for a.a. t ∈ (−τ, T − τ )

Proof. By multiplying the equation by ∂ t z and integrating over Q t , we obtain Qt a|∂ t z| 2 + 1 2 Ω |∇z(t)| 2 = 0 for every t ∈ [0, T ].

Since a is strictly positive, we deduce that both ∂ t z and ∇z vanish a.e. in Q, whence z is a constant. As z(0) = 0, we conclude that z = 0. , we deduce that n ∈ C 0 (Q), whence also κ(n) ∈ C 0 (Q). Now, for p ∈ [2, +∞), we consider the problem of finding w satisfying

e. in Q and w(0) = K(n 0 ) .

(5.13)

Thanks to Proposition 5.1, u is a solution. On the other hand, κ(n) is continuous as just observed, the right-hand side of the equation is bounded and K(n 0 ) belongs to W 2,∞ (Ω) and satisfies ∂ ν K(n 0 ) = 0 on Γ (recall that K is of class C 2 with bounded derivatives by (2.23) and (5.4)). Therefore, we can apply [10, Thm. 2.1] and ensure the existence of a unique w ∈ W 1,p (0, T ; L p (Ω))∩L p (0, T ; W 2,p (Ω)) that satisfies (5.13) and the homogeneous Neumann boundary condition. Since p ≥ 2, (5.12) holds as well and the assumptions of Lemma 5.2 with a = 1/κ(n) are fulfilled by w − u. We conclude that w = u, so that u ∈ W 1,p (0, T ; L p (Ω)) ∩ L p (0, T ; W 2,p (Ω)).

Since p is arbitrary in [2, +∞) , Ω is bounded, n = K −1 (u) and K −1 is of class C 2 with bounded derivatives, we deduce that (2.25) holds for n. Now, we consider the components s, i and h. We observe that each of equations (2.18)-(2.20) and the corresponding boundary and initial conditions have the form ∂ t w − div(κ(n)∇w) = f with ∂ ν w = 0 and w(0) = w 0 (5.14)

where f ∈ L ∞ (Q) and w 0 ∈ W 2,∞ (Ω) with ∂ ν w 0 = 0. On account of the already proved regularity of n, we can rewrite the equation in (5.14) in the non-divergence form, namely ∂ t w − κ(n)∆w + κ ′ (n)∇n · ∇w) = f (5.15) and apply [10, Thm. 2.1] once more with p ≥ 2 (using the 4p-summability of ∇n). Then the homogeneous Neumann problem for (5.15) complemented with the initial condition w(0) = w 0 has a unique solution w ∈ W 1,p (0, T ; L p (Ω)) ∩ L p (0, T ; W 2,p (Ω)). Clearly, such a w belongs to H 1 (0, T ; H) ∩ L 2 (0, T ; W ) and satisfies (5.14). Since (5.14) has a unique solution in H 1 (0, T ; H) ∩ L 2 (0, T ; W ), we deduce that w ∈ W 1,p (0, T ; L p (Ω)) ∩ L p (0, T ; W 2,p (Ω)) whenever w belongs to H 1 (0, T ; H) ∩ L 2 (0, T ; W ) and satisfies (5.14) . This is the case for s, i and h. Since p ≥ 2 is arbitrary and Ω is bounded, (2.25) is completely proved.

Conclusion of the proof of Theorem 2.4. We pick any two solutions (n j , s j , i j , h j ), j = 1, 2, and prove that they are the same. Thanks to Proposition 5.3, they satisfy (2.25) . In particular, the gradients of all the components belong to L 4 (0, T ; L ∞ (Ω)) since W 2,4 (Ω) ⊂ W 1,∞ (Ω). For simplicity, we use the rule of Notation 4.4 concerning the constants: the symbol c stands for possible difference constants that depend on the structure, the data, some norms of the solutions we have fixed and the constants appearing in (5.4) (which are chosen after fixing the solutions) and the consequent Lipschitz constants of the nonlinearities. As in the proof of Theorem 2.3, we set for brevity n := n 1 −n 2 , s := s 1 −s 2 , i := i 1 − i 2 and h := h 1 − h 2 , write equations (2.17)-(2.20) for both solutions and test the differences by n, s, i and h, respectively. As for the first equation, we have While we can simply use the inequality κ(n 1 ) ≥ κ * on the left-hand side, the true novelty is the first term on the right-hand side. We treat it by owing to the Hölder and Young inequalities this way Qt (κ(n 1 ) − κ(n 2 ))∇n 2 · ∇n ≤ c t 0

The other terms on the right-hand side can be dealt with as in the above proof. By combining and rearranging, we easily deduce that Ω |n(t)| 2 ≤ c t 0 ∇n 2 (t ′ ) 2 ∞ + 1 n(t ′ ) 2 + i(t ′ ) 2 dt ′ . Moreover, we can use the inequality κ(n 1 ) ≥ κ * on the left-hand side and rearrange also in this case. By treating equations (2.19)-(2.20) in the same way and summing up, we conclude that Ω |n(t)| 2 + |s(t)| 2 + |i(t)| 2 + |h(t)| 2 ≤ t 0 ψ(t ′ ) n(t ′ ) 2 + s(t ′ ) 2 + i(t ′ ) 2 + h(t ′ ) 2 dt ′ for every t ∈ [0, T ] with a function ψ ∈ L 1 (0, T ). Thus, the Gronwall lemma yields n = s = i = h = 0, and the proof is complete.

Control with uncertain data of socially structured compartmental epidemic models

Nonlinear Differential Equations of Monotone Types in Banach Spaces

A multiscale model of virus pandemic: heterogeneous interactive entities in a globally connected world

Occurrence vs. absence of taxisdriven instabilities in a May-Nowak model for virus infection

Chemotaxis and Cross-diffusion Models in Complex Environments: Models and Analytic Problems Towards a Multiscale Vision

Special Issue on "Mathematics Towards COVID19 and Pandemic

Propagation of epidemics along lines with fast diffusion

Hyperbolic compartmental models for epidemic spread on networks with uncertain data: application to the emergence of COVID-19 in Italy

Continuous-time stochastic processes for the spread of COVID-19 disease simulated via a Monte Carlo approach and comparison with deterministic models

Optimal L p -L q -estimates for parabolic boundary value problems with inhomogeneous data

Degenerate Parabolic Equations

Spread and dynamics of the COVID-19 epidemic in Italy: effects of emergency containment measures

Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy

Adaptive mesh refinement and coarsening for diffusion-reaction epidemiological models

Assessing the spatio-temporal spread of COVID-19 via compartmental models with diffusion in Italy, USA, and Brazil

Delay differential equations for the spatially resolved simulation of epidemics with specific application to COVID-19 Math

Bayesian-based predictions of COVID-19 evolution in Texas using multispecies mixture-theoretic continuum models

Is it safe to lift COVID-19 travel bans? The Newfoundland story

Quelques méthodes de résolution des problèmes aux limites non linéaires

SUIHTER: a new mathematical model for COVID-19. Application to the analysis of the second epidemic outbreak in Italy

Simulating the spread of COVID-19 via a spatiallyresolved susceptible-exposed-infected-recovered-deceased (SEIRD) model with heterogeneous diffusion

Diffusion-reaction compartmental models formulated in a continuum mechanics framework: application to COVID-19, mathematical analysis, and numerical study

System inference for the spatio-temporal evolution of infectious diseases: Michigan in the time of COVID-19

Boundedness in a chemotaxis-May-Nowak model for virus dynamics with mildly saturated chemotactic sensitivity

An agent-based computational framework for simulation of global pandemic and social response on planet X

We treat each term of this equality, separately. Since n = K −1 (u) and the derivative of K −1 is bounded from below by 1/κ * , we haveAs for the second term, we recall that u = ∂ t U, whenceFinally, the Young inequality yieldsBy collecting all this, rearranging and integrating over (−τ, T − τ ) we obtainNow, we observe that the last term can be estimate uniformly with respect to τ by the norm of f and n 0 in L 2 (0, T ; H) and W 2,∞ (Ω), respectively, and that U(0) − U(−τ ) = τ K(n 0 ). Therefore, we conclude that ∂ t u ∈ L 2 (0, T ; H). As n = K −1 (u) and K −1 is Lipschitz continuous, we infer that ∂ t n ∈ L 2 (0, T ; H), so that (5.7) is proved. This conclusions also imply that ∂ t u = κ(n)∂ t n a.e. in Q. Next, we come back to (5.10) written on (0, T ) and recall that f ∈ L ∞ (Q). Thus, (5.8) follows. Finally, we can rewrite (5.10) on (0, T ) this wayThis implies that ∆u ∈ L 2 (0, T ; H) and that 1 κ(n) ∂ t u − ∆u = f a.e. in Q and ∂ ν u = 0 on Γ.Then (5.9) follows. Then z = 0.