Yann Demichel

Anne Estrade

Marie Kratz

Gennady Samorodnitsky

How fast can the chord-length distribution decay?

Keywords: 1991 Mathematics Subject Classification. 60K05, 60D05, 60G10, 60G55, 60G70 chord-length, crossings, Gaussian fields, bi-phasic medium, tail of distribution. 1

The modelling of random bi-phasic, or porous, media has been, and still is, under active investigation by mathematicians, physicists or physicians. In this paper we consider a thresholded random process X as a source of the two phases. The intervals when X is in a given phase, named chords, are the subject of interest. We focus on the study of the tails of the chord-length distribution functions. In the literature concerned with real data, different types of the tail behavior have been reported, among them exponentiallike or power-like decay. We look for the link between the dependence structure of the underlying thresholded process X and the rate of decay of the chord-length distribution. When the process X is a stationary Gaussian process, we relate the latter to the rate at which the covariance function of X decays at large lags. We show that exponential, or nearly exponential, decay of the tail of the distribution of the chord-lengths is very common, perhaps surprisingly so.

Introduction

Studying porous media, such as human bones, food, rocks, etc., leads naturally to D-dimensional Boolean models describing presence or absence of material. Mathematically, a Boolean model is a function f : R D → {0, 1}, when the part of the space where the function f takes value zero represents the "empty" part (lack of material, or "pore"), while the part of the space where the function f takes value one represents the "full" part (presence of material, or "matrix"). A Boolean model is often chosen to be stochastic, and a possible stochastic Boolean model is obtained by thresholding a random field (X t ) t∈R D at a given level γ :

f (t) = 1I (γ,∞) (X t ) = 1 if X t > γ , 0 otherwise.
This procedure is very commonly considered by physicists (see [START_REF] Roberts | Transport properties of heterogeneous materials derived from Gaussian random fields: bounds and simulation[END_REF][START_REF] Nott | Multi-phase image modelling with excursion sets[END_REF][START_REF] Bron | Modelling a food microstructure by random sets[END_REF] for instance).

Here we adopt a one dimensional point of view: we draw test lines through the random medium (X(x), x ∈ R 3 ), and for any line ∆, we identify (X(x), x ∈ ∆) with a process X = (X t , t ∈ R).

The successive intervals with f (t) = 0 (respectively f (t) = 1) are called chords. The chords have been studied in the physics and mathematics literature; see for instance [START_REF] Torquato | Random Heterogeneous Materials[END_REF][START_REF] Roberts | Chord-distribution functions of three-dimensional random media: approximate first-passage times of Gaussian processes[END_REF][START_REF] Stoyan | Stochastic Geometry and its Applications[END_REF] and references therein. In particular, the chord-lengths have been investigated. In a previous paper ( [START_REF] Estrade | Chord-distribution functions and Rice formulae[END_REF]) we defined analytically the chord-length distribution functions. In this paper we focus on their rate of decay.

More precisely, let X = (X t ) t∈ R be a continuously differentiable real strictly stationary process defined on some probability space (Ω, F, P). If the process has a finite variance, we denote its covariance function by ρ. All processes considered in this paper have a finite variance and most of them will be Gaussian.

Assuming the derivative of the process does not vanish on intervals of positive lengths, the chord-lengths are well defined by L 0 = inf{s > 0 : X s = γ} and, for k 0, L k+1 = inf{s > 0 : X s+L 0 +...+L k = γ} ,

where γ is a given level (threshold).

Empirically, both exponential-like and power-like rate of decay of the distribution of the chord-lengths have been observed on real data; see for instance the chapter 2 in [START_REF]Handbook in Porous Solids[END_REF]. Note that this refers to the rate of decay of the probability that a chord is very long. In this paper we investigate the effect of the memory in the thresholded process X on the rate of decay of this probability. In the case of stationary Gaussian processes the memory is measured by the covariance function ρ. Section 1 investigates the tail of the distribution of the chord-lengths statistically, on simulated data. We attempt to discriminate between light or heavy tails (of the chord-lengths) both using the Mean Excess Plot method as a graphical method, and estimating the shape parameter of the associated Generalized Pareto Distribution. The numerical results obtained there motivate a probabilistic analysis developed in the next section. When dealing with stationary Gaussian processes having vanishing memory (ρ → 0 at infinity), as it is the case in the numerical examples, one of our main results (Theorem 2.4, (ii)) shows that for a thresholded Gaussian process with an exponentially fast decreasing covariance function, the tail of the distribution of the chord-lengths is decaying exponentially fast as well. Perhaps, even more surprisingly, we prove that, for underlying Gaussian processes whose covariance function is only assumed to decay to zero at any speed at all, the chord-length distribution decays faster than any negative power function (Theorem 2.4, (i)). This is also true for all r-mixing processes (Theorem 2.6). These theoretical results are proved in Section 2. In order to make them more intuitive, we first study the chordlength distribution decay in the simple case whenever the thresholded process is m-dependent (Proposition 2.3).

Statistical analysis of simulated data

In this section we aim to investigate statistically the behavior of the tail of the distribution of the chord-lengths induced by thresholding a stationary Gaussian process. The purpose here is both to illustrate and motivate the bounds of the next section. To generate samples of chord-lengths, we simulate the underlying Gaussian process. We consider the first two chord-lengths, two different thresholds and three different types of covariance functions of the underlying process. [All Matlab codes, samples and outputs are available on the webpage of Y. Demichel.]

1.1. Simulation.

We start with outlying the simulation procedure. Given the covariance function ρ of a 0 mean unit variance Gaussian process X, we simulate the process on a discrete subset {t 1 , . . . , t n } of a compact interval I. These discrete observations will be used to determine approximately the chords completed within the interval I. Note that we need to simulate a mean 0 and variance 1 Gaussian vector (X(t 1 ), . . . , X(t n )) with a covariance matrix R = (ρ(|t it j |)) 1 i,j n . As usual, the key is to get the square root of the covariance matrix.

Since the Cholesky decomposition method is very expensive, we use the Circulant Embedding Matrix method (see [START_REF] Dietrich | Fast and exact simulation of stationary gaussian processes through circulant embedding of the covariance matrix[END_REF]). Recall that in this approach the covariance matrix R is embedded into a circulant matrix C whose eigenvectors are computed with a Fast Fourier Transform. The square root R 1/2 is built from these eigenvectors. This method works if and only if the minimal eigenvalue v -of C is positive and this property is difficult to ascertain a priori. In practice we choose an interval I and a finite grid {t 1 , . . . , t n } ⊂ I and compute the matrix C and the minimal eigenvalue v -. These choices are crucial since a bad choice may yield an untractable circulant matrix C. In order to have a good sample of chords, we have to consider both a large interval I (to ensure that large chords are not missing) and a fine grid {t 1 , . . . , t n } (to ensure that small chords are not missing). In any case, we will miss all the chords that do not fall within the compact interval I.

Gaussian processes with various covariance functions have been simulated using this method and its extensions; see for instance [START_REF] Gneiting | Fast and Exact Simulation of Large Gaussian Lattice Systems in R 2 : Exploring the Limits[END_REF][START_REF] Schlather | Simulation and Analysis of Random Fields[END_REF][START_REF] Schlather | Reference manual for the Random Fields packages in R[END_REF]. We have chosen three types of covariance functions ρ according to their speed of decay, to see if that would imply different types of chord-length tail behavior. For further illustration, we consider the first two chord-lengths.

The different covariance functions ρ and certain related properties are described in Table 1. We have chosen the time interval I = [0, 3] and a grid of 6000 points. In each case the Circulant Embedding Matrix method works since the minimum eigenvalue is positive, thus allowing us to obtain samples of the first two chord-lengths L 0 and L 1 for a specified threshold. We perform statistical analysis on the samples of the first two chord-lengths for each of the above three covariance functions. We use two different thresholds, γ = 0 and γ = 1. The approach we use is common in the Extreme Value Theory (see for instance [START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF]). First, we use a graphical method, the Mean Excess Plot (MEP), to try to judge whether the chord-lengths have light (i.e. exponentially fast decreasing) or heavy (i.e. hyperbolically fast decreasing) tails. Next, we fit a Generalized Pareto Distribution (GPD) to the upper part of the chordlength distributions. The estimated from the considered sample parameters of the latter will then indicate either light or heavy tails of the chord-length distributions.

Even though we presently apply these two methods to simulated data, the same approach could be used on real data as well.

A graphical method: the Mean-Excess Plot.

Recall the excess cumulative distribution function (cdf) F u of a random variable X over a threshold u ∈ R is defined in the Peak Over Threshold (POT) approach as the cdf of Xu conditioned on X > u, namely

F u (x) = P (X -u x | X > u) , x 0 .
The corresponding mean excess function e of X is defined by e(u) = E(Xu | X > u), whenever it exists.

The plot of the mean excess function e is a useful graphical tool to help distinguishing between heavy and light tails. For example, e(u) = 1/λ for any u if X is exponentially distributed with parameter λ; heavy-tailed distribution functions have a mean excess function tending to infinity, typically along an asymptotically straight line; distribution functions with tails decaying faster than exponentially fast are characterized by a mean excess function tending to 0.

In practice, one uses the empirical mean excess plot

{ (X k,n , e n (X k,n )) : k ∈ {1, . . . , n -1} } ,
where X 1,n • • • X n,n are the order statistics of an n-sample (X i ) 1 i n , and e n (u) is the empirical mean excess function defined using the empirical cumulative distribution function by

1 N u j∈In(u) (X j -u) with I n (u) = {j : 1 j n, X j > u} and N u = Card(I n (u)) .
In the case of u equal to one of the order statistics, this is equivalent to

e n (X k,n ) = 1 n -k n j=k+1 X j,n -X k,n .
Fitting a GPD to the excesses over a threshold. Pickands proved in [START_REF] Pickands | Statistical inference using extreme order statistics[END_REF] that, for sufficiently high threshold u, the excess cdf F u of any random variable X in a domain of attraction of an extreme value distribution can be well approximated by a GPD G ξ,σ(u) , with a shape parameter ξ and scale parameter σ = σ(u) > 0:

G(y) = G ξ,σ(u) (y) =        1 -1 + ξ y σ(u) -1/ξ if ξ = 0, 1 -exp - y σ(u) otherwise,
where y 0 if ξ 0 and 0 y -σ(u)/ξ if ξ < 0. Most of the "textbook" random variables are in the domain of attraction of some extreme value distribution, and so the above approximation of the excess cdf is very general. The shape parameter ξ > 0 arises when X is heavy tailed, and ξ ≤ 0 corresponds to light tails. Therefore, the parameters of the fitted GPD distribution provide information on the tails of X.

An important question is how to select an appropriate high threshold u; we choose it by plotting the empirical mean excess function and choosing u in the range where the latter appears to be linear or stable.

The parameters of a GPD can be estimated via different methods. We will use the method of moments; see [START_REF] Hosking | Parameter and quantile estimation for the Generalized Pareto distribution[END_REF]. If (Y j ) 1 j Nu denote the excesses over a given threshold u in a given sample, then the moments estimators of the parameters ξ and σ(u) of the approximating GPD are given, respectively, by

ξ = 1 2 1 - Y 2 S 2 Y and σ = σ(u) = Y 1 2 + Y 2 S 2 Y , (2)
where Ȳ and S 2 Y are the sample mean and variance of the excesses:

Y = 1 N u Nu i=1 Y i and S 2 Y = 1 N u -1 Nu i=1 Y i -Y 2 .
Provided that the shape parameter satisfies ξ < 1/4, it can be shown by standard methods that the random vector ( σ, ξ) is asymptotically normal with covariance matrix A satisfying, as the sample sizes increases,

N u A ∼ Γ = (1 -ξ) 2 (1 -2ξ)(1 -3ξ)(1 -4ξ) (a ij ) 1 i,j 2 , with a 11 = 2σ 2 (u)(1 -6ξ + 12ξ 2 ), a 22 = (1 -2ξ) 2 (1 -ξ + 6ξ 2 )
and

a 12 = a 21 = σ(u)(1 -2ξ)(1 -4ξ -12ξ 2 ),
from which a confidence interval with asymptotic confidence level α can be deduced:

σ ξ + 1 N u Γ 1/2 q((1 -α)/2) q((1 -α)/2) σ(u) ξ σ ξ + 1 N u Γ 1/2 q((1 + α)/2) q((1 + α)/2) (3)
with q(x) denoting the xth quantile of the standard normal distribution.

1.3. Application to the chord-lengths.

We generated samples of the first chord-lengths (L 0 i ) 1 i n and of the second chord-lengths (L 1 i ) 1 i n of the size n = 10000 each. In each case we started by plotting the empirical mean excess function, in order to judge whether it appears to increase linearly for large levels, or to decay to zero. This was done for each of the three types of covariance and for the two chosen thresholds γ = 0 and γ = 1. Next, on each such mean excess plot, we selected a level u in the range where e n looks approximately linear or stable. It is well known that selecting a proper threshold u is not an easy task as it implies a balance between bias and variance: too high a value of u leads to too few exceedances and, consequently, high variance of the estimators, whereas too small a value of u increases the bias of the estimators. The standard practice is to adopt as low a threshold as possible subject to the limiting GPD model providing a reasonable approximation to the empirical tail. We assess this graphically as well.

Having selected a level u, we estimate the corresponding GPD parameters as in [START_REF] Bradley | Basic properties of strong mixing properties[END_REF], with the associated asymptotic confidence intervals (CI), given in (3), at the confidence level 99%.

The results obtained for the two first chord-lengths, L 0 and L 1 , are summarized in the Table 2 for the threshold γ = 0, and in the Table 3 for the threshold γ = 1. The corresponding empirical mean excess plots and the comparison of the tail of the empirical CDG and the tail of the approximating GPD appear in Figures 1 to 4. The empirical quantile function of order α is denoted by q n (α). It is obvious that the empirical mean excess functions (the left plots of the above figures) appear to be decaying to zero, which is consistent with light tails of the chord-length distributions. Furthermore, the estimated values of the shape parameter ξ of the approximating GDP are all negative. Even though most of the 99% confidence intervals for ξ contain the origin, it is clear that the estimated shape parameter points towards light tails of the chord-length distributions as well. The cdf of the sample and its approximating GPD are represented on the right side of the previous figures. The non-linear behavior of the log-log plots rules out the hypothesis of a power-law decay. The results do not seem to differ significantly for the two chord-lengths, for the two chosen thresholds, and for the three covariance functions ρ. We note that a statistical test of an exponential decay could also be applied to the tails of the chord-length distributions; see e.g. [START_REF] Chaouche | A statistical test procedure for the shape parameter of a generalized Pareto distribution[END_REF].

(b) Chord-length L 1 with γ = 0. Chord-length L 1 Ex. u α s.t. u = q n (α) N u ξ σ CI(ξ) CI(σ) n • 1 0.
(b) Chord-length L 1 with γ = 1. Chord-length L 1 Ex. u α s.t. u = q n (α) N u ξ σ CI(ξ) CI(σ) n • 1

Conclusion:

our empirical analysis appears to indicate light tails of the chord-length, and the result does not seem to be sensitive to the choice of a threshold or a covariance function.

To understand this phenomenon we turn now to a probabilistic analysis.

Theoretical bounds for the tail of chord-length distributions

In this section we derive upper bounds for the tail of the chord-length distribution for certain families of stochastic processes, mostly stationary Gaussian processes, with appropriate assumptions on the memory of the process (on the covariance function if the process is Gaussian). We will see that under mild assumptions a faster than polynomially fast rate of decay of the chord-length distribution is obtained while stronger assumptions lead to an exponentially fast rate of decay. The bounds we obtain provide an explanation to the numerical results discussed in Section 1.

We start with introducing the terminology related to the speed of decay of the tail of distributions.

Definition. Let F be a mapping from [0, ∞) into [0, ∞). We say that F decays exponentially fast if

∃ θ > 0 , ∃ T 0 > 0 , ∀ t T 0 , F (t) e -θt (4)
and decays faster than polynomially fast if

∀ β > 0, lim t→∞ t β F (t) = 0 . (5)
Clearly, the two notions are not mutually exclusive.

Let X = (X t , t 0) be a continuously differentiable strictly stationary process with mean 0 and variance 1.

There is a natural dichotomy between the behavior of chord-lengths if the starting point X 0 is below or above level γ. To account for this, we will write for any chord-length L k and any t 0, P(L k > t) = P(L k > t, X 0 > γ) + P(L k > t, X 0 γ) . Then one estimates each term in the right hand side, essentially in the same way. Moreover, we invariably use a discretization of time.

General inequalities and example.

We start with two general inequalities. The first inequality will allow us to derive asymptotic upper bounds, under the stationary probability measure P, on the tail of the distribution of any chord-length L k from the ones obtained for L 0 . The second inequality will make it possible to extend the asymptotic upper bounds on the tail of the chord-length distribution from the stationary probability P to the Palm probability measures. Finally, we exhibit a simple example in which the underlying process X has memory that does not last longer than m units of time. This example already contains the key idea used in the sequel.

Rate of decay of the tail of distribution of the k-th chord-length.

Proposition 2.1. For any k 0 and any t 0, one has

P(L k > t) 2 k P L 0 > t (k + 1)! .
In particular, if the tail of the distribution of L 0 decays exponentially fast, or faster than polynomially fast, respectively, then so does the tail of the distribution of L k for any k 0.

Proof. The inequality can be easily be proved by induction in k, beginning with k = 0, and using the following decomposition:

P(L k > t) = P(∃ j ∈ {0, . . . , k -1} , L j > t/(k + 1) , L k > t) + P(∀ j ∈ {0, . . . , k -1} , L j t/(k + 1) , L k > t) k-1 j=0 P(L j > t/(k + 1)) + P ∩ s∈(kt/(k+1),t) (X s > γ) ∪ ∩ s∈(kt/(k+1),t) (X s < γ)
and noticing that, by stationarity, the last term in the right hand side above is less or equal to P(L 0 > t/(k + 1)).

Therefore, from now on, we will derive bounds only for the tail of the distribution of the initial chord-length L = L 0 .

Chord-lengths under the Palm probability measures.

In practice it is of special interest to study the chord-length distribution from the moment a level crossing occurs. For that purpose one introduces the so-called Palm probability measures. For instance, to study the chord-length distributions occurring after an upcrossing of the level γ, we introduce the Palm probability measure (see for instance [START_REF] Baccelli | Elements of queueing theory[END_REF][START_REF] Cramér | Stationary and Related Stochastic Processes[END_REF]) defined by

P 0 + (A) = 1 µ lim τ →0 + 1 τ P (A ∩ (U (-τ, 0) 1)) , A ∈ F ,
where U (s, t) denotes the number of upcrossings in the interval (s, t) of the level γ by the process X. Further, µ = E[U (0, 1)] is assumed to be finite; recall (see [START_REF] Cramér | Stationary and Related Stochastic Processes[END_REF]) that under our assumptions,

P(U (0, t) 1) = µt + o(t) as t → 0 . ( 6
)
The Palm probability measures P 0 -, describing the behavior of the process after a downcrossing, and P 0 , describing the behavior of the process after a crossing, are defined analogously.

The distributions of L under P and under P 0 + are linked by

P 0 + (L > t) = - 1 µ ∂ ∂t P(X 0 < γ, L > t) ; (7)
see e.g. [START_REF] Baccelli | Elements of queueing theory[END_REF][START_REF] Cramér | Stationary and Related Stochastic Processes[END_REF][START_REF] Stoyan | Stochastic Geometry and its Applications[END_REF].

Even though the distributions of the chord-lengths under the Palm probability measures are of main interest, they are difficult to evaluate. Fortunately, we have the following:

Proposition 2.2. For any t > 0, one has

P 0 + (L > 2t) 1 tµ P(L > t) .
In particular, if the tail distribution of L under the stationary probability measure P decays exponentially fast, or faster than polynomially fast, respectively, then the same is true under the Palm probability measure P 0 + .

Proof.

Writing P 0 + (L > 2t) 1 t 2t t P 0 + (L > s)
ds and using ( 7) allow to conclude.

A similar argument gives the corresponding bounds for the other Palm probabilities:

P 0 -(L > 2t) 1 tµ P(L > t) and P 0 (L > 2t) 1 2tµ P(L > t) .

2.1.3.

Example: an m-dependent process X. Recall that a stochastic process X is m-dependent if X s and X t are independent whenever |t -s| > m. For second order stationary processes, the m-dependence implies that the covariance function of the process vanishes after lag m, while for stationary Gaussian processes the converse statement is true as well. Processes with such covariance functions are commonly used in simulation (see e.g. [START_REF] Schlather | Reference manual for the Random Fields packages in R[END_REF]). It is the case for our example n • 1 in the simulation part (see §1.1).

Proposition 2.3. Assume that X is an m-dependent stationary process such that P(X 0 > γ) > 0 and P(X 0 < γ) > 0. Then the tail of the distribution of the chord-length L decays exponentially fast.

Proof. Let t > m ′ > m and n = [t/m ′ ] 1. Clearly,

P(L > t, X 0 < γ) P(L > nm ′ , X 0 < γ) n k=0 P(X km ′ < γ) = P(X 0 < γ) n+1 e t log P(X 0 <γ) m ′
.

We deduce that

P(L > t) 2e -θt for all t > m ′ , with θ = 1 m ′ max {| log P(X 0 < γ)|, | log P(X 0 γ)|} .
This implies the exponentially fast rate of decay of the tail.

We mention that the above result remains valid when the m-dependence property of the process X is replaced by certain strong mixing properties such as ψ-mixing (see [START_REF] Bradley | Basic properties of strong mixing properties[END_REF] for a precise definition).

We now consider the case of stationary Gaussian processes under weaker assumptions than the m-dependence considered above.

Gaussian processes with vanishing memory.

Assume that the covariance function ρ of a stationary Gaussian process X tends to 0 at infinity; we will prove that the tail of the distribution of the associated chord-length L decays faster than polynomially fast. When having a fast rate of decay for ρ, the distribution L will decay even faster. The proof is based on the application of Slepian's lemma, introducing a new Gaussian process whose covariance is compared with ρ. Without loss of generality, we may and will assume throughout the section that the stationary Gaussian process has zero mean and unit variance.

Theorem 2.4. Let X be a stationary Gaussian process with covariance function ρ such that ρ(t) → 0 as t → ∞.

(i) The tail of the distribution of the chord-length L decays faster than polynomially fast. (ii) Moreover if ρ decays exponentially fast, then the distribution of L decays exponentially fast as well.

Proof.

(i) Choose 0 < a < 1, and let T be so large that ρ(s) a for all s ≥ T . We use the discretization

P(L > t, X 0 > γ) P X jT > γ, j = 0, 1, . . . , [t/T ] . By Lemma 2.5 below, lim sup t→∞ P(L > t, X 0 > γ) log t lim sup t→∞ P X jT > γ, j = 0, 1, . . . , [t/T ] log[t/T ] log[t/T ] log t - 1 -a a .
Letting a ↓ 0 we obtain lim sup

t→∞ P(L > t, X 0 > γ) log t = -∞ .
Applying the above to the process -X we obtain also lim sup

t→∞ P(L > t, X 0 < γ) log t = -∞ ,
and the two statements together give us the claim of the theorem.

Lemma 2.5. Let 0 < a < 1 and let Y 1 , Y 2 , . . . be a centered unit variance Gaussian process such that Cov(Y i , Y j ) a for all i = j. Then for any γ ∈ R,

lim sup n→∞ log P(Y 1 > γ, . . . , Y n > γ) log n - 1 -a a . (8)
Proof. Let W 0 , W 1 , W 2 , . . . be i.i.d. standard normal random variables, and let

Z j = a 1/2 W 0 + (1 -a) 1/2 W j , j = 1, 2, . . .. Then Z 1 , Z 2 , .
. . is a discrete time centered unit variance Gaussian process such that Cov(Z i , Z j ) = a for all i = j. By the Slepian inequality (see [START_REF] Slepian | The one-sided barrier problem for Gaussian noise[END_REF]) we know that for any n 1 and γ ∈ R, P(Y 1 > γ, . . . , Y n > γ) P(Z 1 > γ, . . . , Z n > γ). Therefore, it is enough to prove [START_REF] Doukhan | Mixing: properties and examples[END_REF] with (Y j ) replaced by (Z j ).

Choose any 0 < θ < 1a a and write

P(Z 1 > γ, . . . , Z n > γ) = ∞ -∞ φ(x) Ψ γ -a 1/2 x (1 -a) 1/2 n dx = x>(2θ log n) 1/2 + x (2θ log n) 1/2 := I 1 (n) + I 2 (n) ,
where φ and Ψ are the density and the tail of a standard normal random variable, respectively. We use the bounds Ψ(x) 1 2 e -x 2 /2 for x ≥ 0, and Ψ(x) > e -(1+ε)x 2 /2 for x large enough, and for any ε > 0. First of all,

I 1 (n) Ψ (2θ log n) 1/2 1 2 e -θ log n = 1 2 n -θ . (9)
On the other hand, selecting 0 < ε < 1 so small that θa(1

+ ε) 1 -a < 1 , (10)
we have for n large enough,

I 2 (n) Ψ γ -a 1/2 (2θ log n) 1/2 (1 -a) 1/2 n = 1 -Ψ a 1/2 (2θ log n) 1/2 -γ (1 -a) 1/2 n 1 -exp - (1 + ε)aθ log n 1 -a n = 1 -n -(1+ε)aθ/(1-a) n ,
and, using [START_REF] Estrade | Chord-distribution functions and Rice formulae[END_REF], we see that for n large enough,

I 2 (n) exp -n 1-(1+ε)aθ/(1-a) /2 = o(n -θ ) .
Combining this bound with (9), we conclude that lim sup

n→∞ log P(Z 1 > γ, . . . , Z n > γ) log n -θ .
Since θ can be taken arbitrarily close to (1a)/a, we obtain (8) for (Z j ) replacing (Y j ), as required.

(ii) Let T > 0 be the positive number from the definition (4) of the exponentially fast decay of ρ, and θ > 0 be the corresponding exponent. We discretize the time parameter of the process X, defining for n 0,

Y n = X(nT ). Then Y = (Y n ) n 0 is a centered unit variance discrete time stationary Gaussian process, whose covariance function ρ Y satisfies ρ Y (k) = ρ(kT ) e -θkT , k = 0, 1, 2, . . . . (11)
For t T , let n = [t/T ] 1, and note that P(L > t) P(Y j > γ , j = 0, . . . , n) + P(Y j < γ , j = 0, . . . , n) .

We will only spell out the procedure for obtaining an upper bound for the first term in the right hand side above. We will prove that there exists α 0 > 0 and N 0 0 such that for all n N 0 , P(Y j > γ , j = 0, . . . , n) e -α 0 n .

This will, clearly, imply that the probability P(L > t) decays exponentially fast.

As previously, we will use the Slepian lemma and introduce another centered unit variance discrete time stationary Gaussian process Z = (Z n ) n 0 , with covariance function ρ Z equal to the upper bound on the covariance function ρ Y in [START_REF] Gneiting | Fast and Exact Simulation of Large Gaussian Lattice Systems in R 2 : Exploring the Limits[END_REF], i.e.

ρ Z (k) = r k 0 , k = 0, 1, 2, . . . , with r 0 = e -θT .

Such process Z does exists; in fact, it can be represented as a causal AR(1) process defined by

Z n+1 = r 0 Z n + ξ n+1 , n 0,
where (ξ n ) n 0 is a Gaussian white noise N (0, 1r 2 0 ). In particular, Z is a Markov process. This property is important since it allows to proceed in a similar way as in the m-dependent case. From the Slepian normal comparison lemma, we know it is enough to prove [START_REF]Handbook in Porous Solids[END_REF] for Z instead of Y .

The threshold γ in ( 12) can be of any sign. Obviously, once we prove the statement for γ < 0, its validity for any other γ will follow. Nonetheless, since there is a particularly simple argument in the case γ > 0 that helps to understand the trick when tackling the case γ 0, we present it first.

If γ > 0, we can use the simple bound

P(Z j > γ , j = 0, . . . , n) P   n j=0 Z j > (n + 1)γ   . (13)
The random variable n j=0 Z j has the normal distribution N (0, σ 2 n ) with σ 2 n n 1 + r 0 1r 0 , so we obtain

P(Z j > γ , j = 0, . . . , n) Ψ nγ 1 -r 0 1 + r 0 1/2
e -α 0 n , for n large enough, using once again the standard upper bound for Ψ.

Note that any α 0 < γ(1r 0 ) 2(1 + r 0 ) can be used above.

When γ 0, the estimate ( 13) is no longer sufficient for our purposes. Let γ ′ = 3 -2γ > 0 and consider the event

A n = {At least ⌈n/3⌉ out of Z 0 , Z 2 , . . . , Z 2n-2 are larger than γ ′ } ,
where ⌈•⌉ denotes the ceiling function.

First, note that P((Z j > γ , j = 0, . . . , 2n -1)

∩ A n ) P   n-1 j=0 Z 2j > (n -⌈n/3⌉)γ + ⌈n/3⌉γ ′   P   n-1 j=0 Z 2j > n   .
This puts us in a situation analogous to the case γ > 0 above, and using the same argument we obtain that for any

α 1 < 1 -r 2 0 2(1 + r 2 0 ) there is N 1 such that for all n ≥ N 1 , P((Z j > γ , j = 0, . . . , 2n -1) ∩ A n ) e -α 1 n . (14)
In order to estimate P (Z j > γ , j = 0, . . . , 2n -1) ∩ A n , let us introduce the sets B n , n ≥ 1 defined by

B n = (b 1 , ..., b n ) ∈ R n : b j > γ, j = 1, .
.., n, and at most ⌈n/3⌉ -1 of the b j 's are larger than γ ′ .

For every n ≥ 1 we define a function on R n by Q n (z 0 , . . . , z 2n-2 ) = P (Z 2j+1 > , j = 0, . . . , n -1 | Z 2j = z 2j , j = 0, . . . , n -1) , (z 0 , . . . , z 2k , . . . , z 2n-2 ) ∈ R n , in the sense of the usual continuous conditional probabilities, and write

P {Z j > γ , j = 0, . . . , 2n -1} ∩ A n = Bn Q n (z 0 , . . . , z 2n-2 )φ n (z 0 , . . . , z 2n-2 )dz 0 . . . dz 2n-2 , (15)
where φ n is the joint pdf of (Z 0 , Z 2 , . . . , Z 2n-2 ).

Given (z 0 , z 1 , . . . , z 2n-2 , z 2n-1 ) such that the vector of the even-numbered coordinates (z 0 , z 2 , . . . , z 2n-2 ) is in B n , the latter vector has at least n -⌈n/3⌉ + 1 coordinates which are strictly smaller than γ ′ . Elementary counting shows that there are at least ⌈n/3⌉ odd numbers j 1 < j 2 < • • • < j ⌈n/3⌉ in the set {1, 3, . . . , 2n -3} such that for every k ∈ {1, . . . , ⌈n/3⌉} z j k -1 < γ ′ and z j k +1 < γ ′ .

It follows from the Markov property of the process Z that

Q n (z 0 , . . . , z 2n-2 ) P Z j 1 > γ, Z j 2 > γ, . . . , Z j ⌈n/3⌉ > γ | Z 2j = z 2j , j = 0, . . . , n -1 = ⌈n/3⌉ k=1 P (Z j k > γ | Z j k -1 = z j k -1 , Z j k +1 = z j k +1 ) .
Recall that all z-values in the conditions are smaller than γ ′ .

Given (Z j

k -1 = z j k -1 , Z j k +1 = z j k +1 ), the random variable Z j k has the normal distribution N (µ k , σ 2 ) with µ k = r 0 1 + r 2 0 (z j k -1 + z j k +1 ) < 2γ ′ r 0 1 + r 2 0 and σ 2 = 1 + r 4 0 (1 + r 2 0 ) 2 > 0 . Therefore Q n (z 0 , . . . , z 2n-2 ) Ψ γ(1 + r 2 0 ) -2γ ′ r 0 1 + r 4 0 ⌈n/3⌉
, and it follows from ( 15) that for all n ≥ 1,

P {Z 0 > γ, . . . , Z 2n-1 > γ} ∩ A n e -α 2 n , with α 2 = Ψ γ(1+r 2 0 )-2γ ′ r 0 √ 1+r 4 0
. Combining this bound with ( 14), shows that there exists N 2 such that for all n ≥ N 2 , P (Z 0 > γ, . . . , Z 2n-1 > γ) e -α 3 n , where we can use any α 3 < min(α 1 , α 2 ). This proves [START_REF]Handbook in Porous Solids[END_REF] with any α 0 < α 3 /2 and some N 0 large enough.

Since a stationary Gaussian process is mixing if and only if its covariance function converges to zero (see [START_REF] Cornfield | Ergodic Theory[END_REF]), Theorem 2.4 (i) states that any mixing stationary Gaussian process has the property that the distribution of L decays faster than polynomially fast. We do not know at the moment if this property holds also for non-Gaussian mixing stationary processes. Nevertheless, the conclusion holds under a stronger dependence condition when assuming a non-Gaussian process is r-mixing (see [START_REF] Doukhan | Mixing: properties and examples[END_REF]).

2.3. r-mixing processes.

Recall that the r-mixing coefficients of a process X are defined by

r T = r(A, B T ) = sup W 1 ∈L 2 X (A) W 2 ∈L 2 X (B T ) |Corr(W 1 , W 2 )| , T > 0 ,
where A = (-∞, 0] and B T = [T, ∞), and for D ⊂ R, L 2 X (D) denotes the closure in L 2 (Ω) of Span{X t , t ∈ D}. See [START_REF] Doukhan | Mixing: properties and examples[END_REF]. In particular, for any T > 0 and events

C 1 , C 2 such that W 1 = 1I C 1 ∈ L 2 X (A) and W 2 = 1I C 2 ∈ L 2 X (B T ), we have r T |P(C 1 ∩ C 2 ) -P(C 1 )P(C 2 )| (P(C 1 )P(C 1 )P(C 2 )P(C 2 )) 1/2 . (16)
We will say that X is an r-mixing process if it satisfies lim

T →+∞ r T = 0 . (17)
Theorem 2.6. Assume that X is a stationary r-mixing process such that P(X 0 > γ) > 0 and P(X 0 < γ) > 0. Then the tail of the distribution of the chord-length L decays faster than polynomially fast.

Proof. Fix T > 0 large enough so that r T < min(P(X 0 < γ), P(X 0 > γ)). Define, for n 0, Y n = X(nT ) and consider the discrete time process Y = (Y n ) n 0 .

For n 0 we apply [START_REF] Roberts | Transport properties of heterogeneous materials derived from Gaussian random fields: bounds and simulation[END_REF] with

C 1 = X(-(2 n -1)T ) > γ, . . . , X(0) > γ and C 2 = X(T ) > γ, . . . , X(2 n T ) > γ . By stationarity, P(Y 0 > γ, . . . , Y 2 n+1 -1 > γ) = P({Y 0 > γ, . . . , Y 2 n -1 > γ} ∩ {Y 2 n > γ, . . . , Y 2 n+1 -1 > γ})
(P(Y 0 > γ, . . . , Y 2 n -1 > γ)) 2 + r T P(Y 0 > γ, . . . , Y 2 n -1 > γ) .

Denoting p n = P(Y 0 > γ, . . . , Y 2 n -1 > γ) for n 0, we see that p n+1 p 2 n + r T p n , n = 0, 1, . . . .

which implies that p n → 0 as n → ∞ (if l = lim n→∞ p n = 0, then l > p 0 which would contradict (p n ) decreasing). Therefore, there exists N 0 0 such that for all n N 0 one has p n r T . Thus p n+1 2r T p n for all n N 0 and, hence, p n c 0 (2r T ) n with c 0 = p N 0 (2r T ) -N 0 > 0 for all n N 0 . Now let t > 2 N 0 T be so large that n = [log 2 (t/T )] N 0 and remark that r T < 1/2. We have, log 2 denoting the logarithm in base 2, P(L > t , X 0 > γ) P(Y 0 > γ, . . . , Y 2 n -1 > γ) c 0 (2r T ) n c ′ 0 t -|log 2 (2r T )| for a certain constant c ′ 0 > 0 depending only on N 0 and T . We deduce that lim t→∞ log P(L > t , X 0 > γ) log t -|log 2 (2r T )| .

Since r T → 0 as T → ∞, we conclude that lim t→∞ log P(L > t , X 0 > γ) log t = -∞ .

Hence t β P(L > t , X 0 > γ) → 0 as t → ∞, for any β > 0. An analogous result holds for P(L > t , X 0 γ). Therefore the distribution of L decays faster than polynomially fast.

Recall that a stationary Gaussian process with a covariance function which decays in a polynomial way is r-mixing (see [START_REF] Cornfield | Ergodic Theory[END_REF]). So an immediate corollary of Theorem 2.6 is the fact that for Gaussian processes with polynomially decaying covariance function, the distribution of the chord-length L decays faster than polynomially fast, a result obtained in Theorem 2.4 (i) under a weaker assumption.

Conclusion

We have established theoretical results on the rate of decay of the tail of the distribution of the chord-lengths, depending on the memory in the thresholded process. In the case of stationary Gaussian processes, the memory is expressed via the covariance function. The results agree with the empirical results obtained by statistical analysis of simulated processes. It shows that, as soon as one deals with a thresholded stationary process whose covariance decays to 0, a rapidly decreasing decay of the chord-length distribution has to be expected. Consequently, it seems hopeless to use such a thresholded process as a model for real data when, for instance, a power-like decay of the chord-length distribution is observed.

We focused on the asymptotic behavior of the tails of the distributions and did not try to derive the best constants in the upper bounds for such tails. It may be possible to refine these bounds and relate them to the threshold. It may further be possible to explore how the bounds change from one chord-length to the next. This will be the subject of a future work. We also plan to investigate what memory properties of the thresholded process imply lower bounds on the tails of the chord-length distributions.

We hope that our work will also contribute to study of the D-dimensional Boolean models for D > 1. For instance, if L * is the spherical contact distance of the phase containing the origin, upper bounds for P(L * > t) can be obtained from our results since L * = sup{R > 0 : B(0, R) ⊂ 0-phase} = inf α∈[0,2π)

L(α) ,

where the 0-phase is the phase containing the origin and L(α) denotes the initial chord-length in the direction α.

10 -Figure 1 .

 101 Figure 1. Length of the chord L 0 at threshold γ = 0 for the three examples of covariance function ρ (from top to bottom: n • 1, n • 2, n • 3). On the left: the empirical mean excess plots; on the right: log-log plots of the tails of the distribution functions (solid line: sample, dashed line: approximating GPD model). The parameters are given in Table 2(a).

Figure 2 .

 2 Figure 2. Length of the chord L 1 at threshold γ = 0 for the three examples of covariance function ρ (from top to bottom: n • 1, n • 2, n • 3). On the left: the empirical mean excess plots; on the right: log-log plots of the tails of the distribution functions (solid line: sample, dashed line: approximating GPD model). The parameters are given in Table 2(b).

Figure 3 .

 3 Figure 3. Length of the chord L 0 at threshold γ = 1 for the three examples of covariance function ρ (from top to bottom: n • 1, n • 2, n • 3). On the left: the empirical mean excess plots; on the right: log-log plots of the tails of the distribution functions (solid line: sample, dashed line: approximating GPD model). The parameters are given in Table 3(a).

Figure 4 .

 4 Figure 4. Length of the chord L 1 at threshold γ = 1 for the three examples of covariance function ρ (from top to bottom: n • 1, n • 2, n • 3). On the left: the empirical mean excess plots; on the right: log-log plots of the tails of the distribution functions (solid line: sample, dashed line: approximating GPD model). The parameters are given in Table 3(b).

Table 1 .

 1 Covariance functions used for simulation.

	Example	Speed of decay	Expression ρ(x) Eigenvalue v -
	n • 1	compact support P 11 (5x)1I [0,1] (5x) 6.8426.10 -10
	n • 2	very fast	exp(-(5x) 2 )	1.1827.10 -11
	n • 3	polynomial	(1 + (6x) 2 ) -4	2.7569.10 -11
	P 11 (x) = 1 -22 3 x 2 + 33 x 4 -77 2 x 5 + 33 2 x 7 -11 2 x 9 + 5 6 x 11 .
	1.2. Statistical analysis.

Table 2 .

 2 Statistical results for L 0 and L 1 with γ = 0. Chord-length L 0 with γ = 0.

				Chord-length L 0
	Ex.	u	α s.t. u = q n (α)	N u	ξ σ	CI(ξ) CI(σ)
	n • 1 0.41	90%	910	-0.0955 0.1816	(-0.2531, 0.0621) (0.1047, 0.2585)
	n • 2 0.90	86%	1374	-0.1530 0.4575	(-0.3301, 0.0241) (0.3169, 0.5981)
	n • 3 0.49	91%	881	-0.1058 0.2286	(-0.2767, 0.0651) (0.1343, 0.3229)

(a)

Table 3 .

 3 Statistical results for L 0 and L 1 with γ = 1.

	35	83%	1696	-0.0005 0.1721	(-0.0953, 0.0843) (0.1302, 0.2140)
	n • 2 1.27	92.5%	746	-0.1911 0.4207	(-0.4526, 0.0704) (0.2238, 0.6176)
	n • 3 0.46	86%	1380	-0.0063 0.2056	(-0.1088, 0.0962) (0.1530, 0.2582)

(a) Chord-length L 0 with γ = 1.

MODAL'X, EA 3454, University Paris Ouest Nanterre La Défense,

avenue de la République, 92001 Nanterre cedex, France.E-mail address: yann.demichel@u-paris10.fr

ESSEC Business School Paris, avenue Bernard Hirsch BP 50105, 95021 Cergy Pontoise Cedex, France. E-mail address: kratz@essec.fr

School of Operations Research and Information Engineering, Cornell University, Ithaca, NY 14853, U.S.A.E-mail address: gennady@orie.cornell.edu

This work was partially supported by the French grant "mipomodim" ANR-05-BLAN-0017. Samorodnitsky's research was also partially supported by the ARO grant W911NF-07-1-0078 at Cornell University.

During the elaboration of this work Yann Demichel was member of MAP5, Univ. Paris Descartes 2 MAP5, UMR CNRS 8145, University Paris Descartes, 45 rue des Saints-Pères, 75270 Paris 06, France.

E-mail address: anne.estrade@parisdescartes.fr


Mauricio Labadie
email: mauricio.labadie@gmail.com

Anna Marciniak-Czochra

A reaction-diffusion model for viral infection and immune response

Keywords: Reaction-diffusion systems, global solutions, uniqueness, asymptotic behavior, virus, interferon, immune response

In this work we extend the ODE model for virus infection and immune response proposed by P. Getto et al (Modelling and analysis of dynamics of viral infection of cells and of interferon resistance, J. Math. Anal. Appl., No. 344, 2008, pp. 821-850) to account for the spatial effects of the processes, such as diffusion transport of virions, biomolecules and cells. This leads to two different nonlinear PDE models, a first one where the cells and the biomolecules diffuse (which we call the reactiondiffusion model) and a second one where only the biomolecules can diffuse (the hybrid model). We show that both the reaction-diffusion and the hybrid models are well-posed problems, i.e., they have global unique solutions which are non-negative, bounded, and depend continuously on the initial data. Moreover, we prove that there exists a "continuous" link between these two models, i.e., if the diffusion coefficient of the cells tends to zero then the solution of the reaction-diffusion model converges to the solution of the hybrid model. We also prove that the solutions are uniformly bounded and integrable for all times. We characterize the asymptotic behavior of the solutions of the hybrid system and present several relations concerning the survivability of viruses and cells. Finally, we show that the solutions of the hybrid model converge to the steady state solutions, which implies that the latter are globally stable. We finish with several numerical simulations performed in Matlab.

1 Introduction

Spatial effects of viral infection and immunity response

When a virion, i.e., an individual viral particle, enters a healthy cell, it modifies the genetic structure of its host. After infection, the altered biochemical machinery of the host starts to create new virions. The virions are then released from the host cell and may infect other cells. However, the infected cell activates intrinsic host defenses, which include, among others, activation of the innate immunity system and release of biomolecules called interferons (IFN), which communicate with the other cells and induce them to deploy protective defenses. The dynamics of such complex virus-host system results from the intra-and extracellular interactions between invading virus particles and cells producing substances, which confer resistance to virus.

These key processes have been addressed by the mathematical model proposed by Getto, Kimmel and Marciniak-Czochra in [START_REF] Getto | Modelling and analysis of dynamics of viral infection of cells and of interferon resistance[END_REF]. The model was motivated by the experiments involving vesicular stomatitis virus (VSV) [START_REF] Duca | Quantifying viral propagation in vitro: Toward a method for characterisation of complex phenotypes[END_REF][START_REF] Lam | Arrested Spread of Vesicular Stomatitis Virus Infections In Vitro Depends on Interferon-Mediated Antiviral Activity[END_REF] and respiratory syncytial virus (RSV). The work of P. Getto et al [START_REF] Getto | Modelling and analysis of dynamics of viral infection of cells and of interferon resistance[END_REF] is focused on the study of the role of heterogeneity of intracellular processes, reflected by a structure variable, in the dynamics of the system and stability of stationary states. It is shown that indeed the heterogeneity of the dynamics of cells in respect to the age of the individual cell infection may lead to the significant changes in the behavior of the model solutions, exhibiting either stabilizing or destabilizing effects.

Another interesting aspect of the dynamics of the spread of viral infection and development of resistance is related to the spatial structure of the system and the effects of spatial processes, such as random dispersal of virions and interferon particles. In a series of experiments on the vesicular stomatitis virus infection [START_REF] Duca | Quantifying viral propagation in vitro: Toward a method for characterisation of complex phenotypes[END_REF][START_REF] Lam | Arrested Spread of Vesicular Stomatitis Virus Infections In Vitro Depends on Interferon-Mediated Antiviral Activity[END_REF], it was observed that the spatial structure of the system may influence the dynamics of the whole cell population. The role of spatial dimension and diffusion transport of virions and interferon molecules were experimentally studied using two type of experiments: a one-step growth experiment in which all cells were infected simultaneously, and a focal-infection experiment in which cell population was infected by a point source of virus. A spreading cicular wave of infection followed by a wave of dead cells was observed. The experiments were performed on two different cell cultures: DBT (murine delayed brain tumor) cells, which respond to IFN and can be activated to resist the replication of viruses, and BHK (baby hamster kidney) cells, which are not known to produce or respond to IFN. In case of focal infection both in DBT and BHK populations spread of infection (rings) was observed. The size of the rings was dependent on the type of the virus (N1, N2, N3, N4 -gene ectopic strains as well as M51R mutant and XK3.1). However, for all virus types, it was observed that in DBT cells the speed of the infection propagation was decreasing with time, while in case of BHK the radius of the infected area was growing linearly in time. Results of the experiments showed that the rate of infectious progeny production in one-step growth experiments was a key determinant of the rate of focal spread under the absence of IFN production. Interestingly, the correlation between one-step growth and focal growth did not apply for VSV strains XK3.1 and M51R in the cells producing IFN. Focal infection in DBT cells led only to the limited infection, the spread of which stopped after a while.

These experiments indicated suitability of focal infections for revealing aspects of virus-cell interactions, which are not reflected in one-step growth curves. Motivated by Duca et al experiments, we devise a model of spatio-temporal dynamics of viral infection and interferon production, which involves virions, uninfected, infected and resistant cells, as well as the interferon. We assume that interferon is produced by infected cells and spread by diffusion to neighboring uninfected cells, making them resistant. At the same time, the virus is spread, also by diffusion and the final outcome is the result of competition beween these two processes.

The original model

The point of departure for this work is a system of nonlinear ordinary differential equations developed by P. Getto et al [START_REF] Getto | Modelling and analysis of dynamics of viral infection of cells and of interferon resistance[END_REF] to model wild-type cells interacting with virions and the consequent interferon-based immune response. The model describes time dynamics of five different types of entities interacting in a culture: three different types of cells (wild-type, infected and resistant) and two different biomolecules (virions and interferons). P. Getto et al [START_REF] Getto | Modelling and analysis of dynamics of viral infection of cells and of interferon resistance[END_REF] considered a culture of wild-type cells infected by virions. When a wild-type cell interacts with a virion it becomes an infected cell. Infected cells produce further virions, but they also release interferons. If an interferon reaches a wild-type cell before a virion does, then this cell becomes a resistant cell. Virions, interferons and infected cells are supposed to have an exponential death rate, whilst there is no death rate associated to wild-type and resistant cells, because they are considered to live longer than infected cells. In the other words, the model describes the time scale of in vitro experiments, which is short comparing to the life span of healthy cells.

Under these hypotheses the following nonlinear ODE model was proposed:

           W = -iW -vW , ← wild-type cells I = -µ I I + vW , ← infected cells R = iW , ← resistant cells v = -µ v v + α v I -α 4 vW , ← virions i = -µ i i + α i I -α 3 iW , ← interferons (1)
where all coefficients are non-negative constants.

Biological hypotheses of the new model

In this work we will modify the ODE model (1) by introducing spatial random dispersion of cells, virions and IFN molecules. We consider two models: In the first model (which we call the reaction-diffusion model) we assume that all cells, virions and interferons diffuse, whilst the second model (which we call the hybrid model) is based on the hypothesis that only virions and interferons diffuse. The diffusion terms are supposed to follow Fick's Law with constant diffusion coefficients and it is modeled by adding Laplacian operators to the ODE system.

Concerning the boundary conditions, we assume that the whole system is isolated within a bounded domain Ω ⊂ R N (N = 2, 3). This implies no-flux boundary conditions, i.e., homogeneous Neumann conditions, on Γ = ∂Ω. It is worth to remark, however, that the results presented here are also valid for other boundary conditions such as homogeneous Dirichlet and Robin (mixed).

Since the cells are far bigger than the interferon molecules and the virions, we can suppose that d is much smaller than both d i and d v . In order to compare the two latter diffusion coefficients, one can recall that interferons are biolomecules and virions in general have several proteins (including DNA or RNA). Under the light of this argument we assume that d v is smaller than d i , but both are of the same order. This leads to the following conditions,

0 < d d v ≤ d i .

The reaction-diffusion (RD) model

Let Ω ⊂ R N (N = 2, 3) be a bounded domain with Lipschitz boundary Γ, and let T > 0. We consider in Ω × [0, T ] a reaction-diffusion (RD) system

           ∂ t W = d∆W -iW -vW , ∂ t I = d∆I -µ I I + vW , ∂ t R = d∆R + iW , ∂ t v = d v ∆v -µ v v + α v I -α 4 vW , ∂ t i = d i ∆i -µ i i + α i I -α 3 iW , (2)
with boundary conditions

           ∇W • n(σ, t) = 0 on Γ × [0, T ], ∇I • n(σ, t) = 0 on Γ × [0, T ], ∇R • n(σ, t) = 0 on Γ × [0, T ], ∇v • n(σ, t) = 0 on Γ × [0, T ], ∇i • n(σ, t) = 0 on Γ × [0, T ], (3)
and initial conditions

           W (x, 0) = W 0 (x) , I(x, 0) = I 0 (x) , R(x, 0) = R 0 (x) , v(x, 0) = v 0 (x) . i(x, 0) = i 0 (x) . (4)

The hybrid model

Since d is much smaller than both d i and d v it is plausible to consider that d = 0. Under this assumption we obtain a hybrid model consisting of PDE equations for the interferons i and virions v and ODE equations for the three types of cells W, I, R. The system takes the form

           ∂ t W = -iW -vW , ∂ t I = -µ I I + vW , ∂ t R = iW , ∂ t v = d v ∆v -µ v v + α v I -α 4 vW , ∂ t i = d i ∆i -µ i i + α i I -α 3 iW , (5)
with boundary conditions ∇v

• n(σ, t) = 0 on Γ × [0, T ], ∇i • n(σ, t) = 0 on Γ × [0, T ]. ( 6
)

Main results

In this Section we formulate main results of the work. The proofs are presented in the following sections.

Existence and uniqueness results

Throughout this work we denote • the norm in L 2 (Ω) and • ∞ the norm in L ∞ (Ω).

We will prove that both the RD and the hybrid models are well-posed problems, i.e., they have unique solutions which are non-negative, uniformly bounded, and depend continuously on the initial data.

Theorem 1 Fix any T > 0. If the initial conditions (4) are non-negative a.e. and lie in L ∞ (Ω) then the RD system (2)-(3) has unique solutions W (x, t), I(x, t), R(x, t), v(x, t) and i(x, t) on Ω × [0, T ]. Moreover, these solutions are non-negative, uniformly bounded, and depend continuously on the initial data.

Theorem 2 Fix any T > 0. If the initial conditions (4) are non-negative a.e. and lie in L ∞ (Ω) then the hybrid system ( 5)-( 6) has unique solutions W (x, t), I(x, t), R(x, t), v(x, t) and i(x, t) on Ω × [0, T ]. Moreover, these solutions are non-negative, bounded, and depend continuously on the initial data.

We also prove that there is a "continuous link" between these two models, as the next result shows. 3), which converges to the solution (W, I, R, v, i) of the hybrid system ( 5)-( 6), in the following sense:

Theorem 3 If d → 0 then there is a subsequence (W d , I d , R d , v d , i d ) of solutions of the RD system (2)-(
• Strongly in L 2 0, T ; L 2 (Ω) . • Weakly in L 2 0, T ; H 1 (Ω) . • Weakly-in L ∞ (Ω × (0, T )).
Moreover, any convergent subsequence (W d , I d , R d , v d , i d ) has the same limit (W, I, R, v, i).

Asymptotic results for the RD system

Theorem 4

1. The solutions W, I, R, v, i of the RD system (2)-( 3) are globally-defined and belong to L ∞ (Ω × (0, ∞)).

2. If W, I, R, v, i are non-negative, steady-state solutions of the RD system (2)-(3) then

W (x) = W 0 ≥ 0 constant, I(x) ≡ 0, R(x) = R 0 ≥ 0 constant, v(x) ≡ 0, i(x) ≡ 0.

Asymptotic results for the hybrid system

Regularity of solutions W, I, R, v, i follows from a classical theory for evolution systems and depends on the regularity of initial conditions (see eg. [START_REF] Rothe | Global Solutions of Reaction-Diffusion Systems[END_REF], for reaction-diffusion systems coupled with ODEs). In the reminder of this paper we assume that the solutions are (at least) C 1 .

Theorem 5 If W, I, R, v, i are non-negative, steady-state solutions of the hybrid system ( 5)-( 6) then

I(x) ≡ 0, v(x) ≡ 0, i(x) ≡ 0.
Moreover, suppose that the initial conditions belong to L ∞ (Ω). Then the solutions of the hybrid system ( 5)-( 6) are globally-defined and have the following asymptotic properties:

1. I(x, t) belongs to L 1 (0, ∞; L 2 (Ω)), i.e., lim t→∞ t 0 Ω I 2 (x, s) dΩ ds < ∞. 2. v(x, t) and i(x, t) belong to L 1 (0, ∞; H 1 (Ω)), i.e., lim t→∞ t 0 Ω v 2 (x, s) dΩ + Ω |∇v(x, s)| 2 dΩ ds < ∞ , lim t→∞ t 0 Ω i 2 (x, s) dΩ + Ω |∇i(x, s)| 2 dΩ ds < ∞ . 3. v(x, t)W (x, t) and i(x, t)W (x, t) belong to L 1 (0, ∞; L 1 (Ω)), i.e., lim t→∞ t 0 Ω v(x, s)W (x, s) dΩ ds < ∞ and lim t→∞ t 0 Ω i(x, s)W (x, s) dΩ ds < ∞ . 4. For any x ∈ Ω, lim t→∞ I(x, t) = 0 , lim t→∞ v(x, t) = 0 , lim t→∞ i(x, t) = 0 . 5. For any x ∈ Ω, W 0 (x) > 0 if and only if lim t→∞ W (x, t) > 0.
Theorem 6 Consider the hybrid system ( 5)-( 6) and suppose that µ v = 0. Then:

1. If v ∞ (x) is a steady-state solution then ∇v ∞ = 0. 2. Define v ∞ (x) := lim sup t→∞ v(x, t). If α v ≥ α 4 µ I then Ω v ∞ (x) dΩ ≥ Ω v 0 (x) dΩ.
In particular, if v 0 ≡ 0 then v ∞ ≡ 0.

Numerical simulations

In order to solve numerically the considered systems we used the method of lines, under which the system of nonlinear partial differential equations was converted to a large system of ordinary differential equations by discretising of Laplacian (three coordinates for the 1D Laplacian and five coordinates in the 2D case). The discretised system of ordinary differential equations was numerically solved using the CVODE package and numerical estimates of the Jacobian matrix. This program offers an implicit method for time discretization, originally developed for stiff problems of ODEs. Space discretization is a gridpoint on a unit interval. The size of the spatial grid is adjusted according to the value of the diffusion coefficient. Time discretization is performed implicitly. Homogeneous Neumann boundary conditions (zero flux) are implemented as a reflection at the boundary (see e.g. [START_REF] Ames | Numerical Methods for Partial Differential Equations[END_REF]). The graphical visualization of the numerical solutions in space and time is realized using Matlab.

The domain is the 2D square [0, 1] × [0, 1]. The initial concentration in the 2D simulations are W 0 = 1 and v 0 = 0.5 on the sub-square [0.4, 0.6] × [0.4, 0.6] and zero otherwise. These initial conditions are exactly those in the numerical simulation of the ODE model in P. Getto et al [START_REF] Getto | Modelling and analysis of dynamics of viral infection of cells and of interferon resistance[END_REF]. For the simulation of the hybrid model we used a 30 × 30 grid in space and 50 time steps. The pictures correspond to the final stage t = 50. We fixed the parameters α 3 = 1, α 4 = 4, µ I = 0.3, µ i = 0.4 and α i = 0.8.

From the numerical simulations we could assess the effect of the spacial structure (i.e. the diffusion) in the virus proliferation:

• When µ v > 0 the viral diffusion plays against virions and helps the wild-type cells. In pictures 3.1-3.3 we can see that, as the diffusion coefficient grows, the concentration of wild-type cells seems to grow.

• When µ v = 0 the viral diffusion has a positive effect on the viral concentration and plays against the wild-type cells. In pictures 3.4-3.6 we can see that, as the diffusion coefficient grows, the concentration of virions grows as well.

As we can see in both cases, the spacial structure (represented by the diffusion term) has crucial effect on the final concentration of wild-type cells and virions. Indeed, if µ v > 0 the diffusion helps the cells and punishes the virions whilst if µ v = 0 it has the opposite effect. [START_REF] Hamdache | On a reaction-diffusion model for calcium dynamics in dendritic spines[END_REF] The fixed point operator and a priori estimates 4.1 Construction of the fixed point operator R Our approach follows the ideas used by K. Hamdache and M. Labadie in [START_REF] Hamdache | On a reaction-diffusion model for calcium dynamics in dendritic spines[END_REF]. Fix (v , i ) ∈ K and set

Define K := (v , i ) ∈ L 2 0, T ; L 2 (Ω) 2 ∩ [L ∞ (Ω × [0, T ])] 2 : v (x, t) ≥ 0, i (x, t) ≥ 0 a.e. in Ω × [0, T ] .
∂ t W -d∆W + i W + v W = 0 in Ω × (0, T ), ∂ t I -d∆I + µ I I -v W = 0 in Ω × (0, T ), ∂ t R -d∆R -i W = 0 in Ω × (0, T ), ( 7
)
∇W • n(σ, t) = 0 on Γ × [0, T ], ∇I • n(σ, t) = 0 on Γ × [0, T ], ∇R • n(σ, t) = 0 on Γ × [0, T ].
For any finite time interval [0, T ], the linear system (7) has a unique solution (W (x, t), I(x, t), R(x, t)), which is non-negative and bounded.

With these functions W (x, t), I(x, t) and R(x, t) set

∂ t v -d v ∆v + µ v v -α v I + α 4 vW = 0 in Ω × (0, T ), ∂ t i -d i ∆i + µ i i -α i I + α 3 iW = 0 in Ω × (0, T ), ( 8
) ∇v • n(σ, t) = 0 on Γ × [0, T ]. ∇i • n(σ, t) = 0 on Γ × [0, T ],
Again, for any finite time interval [0, T ] the linear, uncoupled system (8) has a unique solution (v(x, t), i(x, t)), which is non-negative and bounded.

Our goal is to show that the operator R[(v , i )] := (v, i), defined for the chain of maps (v , i ) → (W, I, R) → (i, v) constructed above, has a fixed point.

Positivity of solutions

From now on we will assume that the coefficients

d, d v , d i , µ I , µ v , µ i , α 3 , α 4
are all positive, and that the initial conditions (4) are non-negative and bounded.

Lemma 1 Let (v , i ) → (W, I, R) → (v, i) be solutions of ( 7)-( 8).

1. If v and i are non-negative and bounded then W , I and R are non-negative and bounded.

2. If W , I and R are non-negative and bounded then v and i are non-negative and bounded.

Proof:

1. The equation for W is

∂ t W -d∆W + (i + v )W = 0 in Ω × (0, T ), ∇W • n(σ, t) = 0 on Γ × [0, T ], W (x, 0) = W 0 (x) ≥ 0 in Ω.
Applying Maximum Principle we obtain that W (x, t) ≥ 0 for all (x, t) ∈ Ω × (0, T ).

Define

Z 1 := W -γ 1 , where γ 1 ∈ R. Then Z 1 solves ∂ t Z 1 -d∆Z 1 + (i + v )Z 1 = -(i + v )γ 1 in Ω × (0, T ), ∇Z 1 • n(σ, t) = 0 on Γ × [0, T ], Z 1 (x, 0) = W 0 (x) -γ 1 in Ω.
Choosing γ 1 = W 0 ∞ and using non-negativity of i and v we obtain

∂ t Z 1 -d∆Z 1 + (i + v )Z 1 ≤ 0 in Ω × (0, T ), ∇Z 1 • n(σ, t) = 0 on Γ × [0, T ], Z 1 (x, 0) ≤ 0 in Ω.
Hence, Maximum Principle implies that Z 1 (x, t) ≤ 0 for all (x, t) ∈ Ω × (0, T ), and in consequence W (x, t) ≤ W 0 ∞ for all (x, t) ∈ Ω × (0, T ).

For I, notice that v and W are non-negative. Therefore, similarly as above Maximum Principle yields that I(x, t) ≥ 0 for all (x, t) ∈ Ω × (0, T ). Now, due to the boundedness of I, the function

Z 2 := I -γ 2 solves ∂ t Z 2 -d∆I + µ I Z 2 = v W -µ I γ 2 in Ω × (0, T ), ∇Z 2 • n(σ, t) = 0 on Γ × [0, T ], Z 2 (x, 0) = I 0 (x) -γ 2 in Ω.
Choosing

γ 2 = max I 0 ∞ , v ∞ W 0 ∞
µ I and applying Maximum Principle, we obtain that W (x, t) ≤ γ 2 for all (x, t) ∈ Ω × (0, T ).

Finally, Maximum Principle implies that R(x, t) ≥ 0 for all (x, t) ∈ Ω × (0, T ). For the boundedness of R, define Z 3 := R -γ 3 (t). It yields

∂ t Z 3 -d∆I = i W -γ 3 (t)
in Ω × (0, T ),

∇Z 3 • n(σ, t) = 0 on Γ × [0, T ], Z 3 (x, 0) = R 0 (x) -γ 3 (0) in Ω. Choosing γ 3 (t) = R 0 ∞ + t i ∞ W 0 ∞
and applying Maximum Principle we obtain that R(x, t) ≤ γ 3 (t) for all (x, t) ∈ Ω × (0, T ).

2. Using the same argument as before we can prove that 0 ≤ v(x, t) ≤ γ 4 and 0 ≤ i(x, t) ≤ γ 5 for all (x, t) ∈ Ω × (0, T ), where

γ 4 = max v 0 ∞ , α v I 0 ∞ µ v , γ 5 = max i 0 ∞ , α i I 0 ∞ µ i .

A priori estimates

Lemma 2 Let (v , i ) → (W, I, R) → (v, i) be solutions of ( 7)- [START_REF] Tartar | An introduction to Navier-Stokes equation and Oceanography[END_REF]. Then the functions W, I, R, v, i belong to L 2 0, T ; H 1 (Ω) .

Proof:

• Multiply ( 7) by W and integrate by parts to obtain 1 2

d dt W 2 + d ∇W 2 = - Ω (i + v )W 2 dΩ .
Since i and v are non-negative, it follows that 1 2

d dt W 2 + d ∇W 2 ≤ 0 .
Integrating over [0, t] yields

W (t) 2 + 2d t 0 ∇W (s) 2 ds ≤ W (0) 2 .
• Multiply ( 7) by I and integrate by parts to obtain 1 2

d dt I 2 + d ∇I 2 + µ I I 2 = Ω v W I dΩ .
Recall the identity

Ω |v W I| dΩ ≤ 1 4ε Ω |v W | 2 dΩ + ε Ω |I| 2 dΩ .
Choosing ε = µ I /2 and using the uniform bounds in Lemma 1, it follows that there is a constant C > 0, depending on the L ∞ norm of the initial data, such that

d dt I 2 + 2d ∇I 2 + µ I I 2 ≤ C .
Here we can keep or discard the term µ I I 2 , leading to two different estimates:

-If we discard the term, the integration over [0, T ] yields

I(t) 2 + 2d t 0 ∇I(s) 2 ds ≤ I(0) 2 + Ct , which implies that I ∈ L 2 0, T ; H 1 (Ω) .
-If we keep the term, after multiplying by e µ I t and integrating over [0, T ],

e µ I t I(t) 2 + 2d t 0 e µ I s ∇I(s) 2 ds ≤ I(0) 2 + C(e µ I t -1) .
Therefore,

I(t) 2 + 2d t 0 e -µ I (t-s) ∇I(s) 2 ds ≤ I(0) 2 e -µ I t + C(1 -e -µ I t ) .
This estimate will be useful for the study of the asymptotic behavior when t → ∞.

• Multiply (7) by R and integrate by parts to obtain 1 2

d dt R 2 + d ∇R 2 = Ω iW R dΩ .
Therefore,

d dt R 2 + 2d ∇R 2 ≤ C + R 2 .
Multiplying by e -t and integrating over [0, t] leads to

R(t) 2 + 2d
t 0 e t-s ∇R(s) 2 ds ≤ R(0) 2 e t + C(e t -1) .

• Repeating the argument for v and i yields

v(t) 2 + 2d v t 0 ∇v(s) 2 ds ≤ v(0) 2 + Ct, v(t) 2 + 2d v t 0 e -µv(t-s) ∇v(s) 2 ds ≤ v(0) 2 e -µvt + C(1 -e -µvt ), i(t) 2 + 2d i t 0 ∇i(s) 2 ds ≤ i(0) 2 + Ct, i(t) 2 + 2d i t 0
e -µi(t-s) ∇i(s) 2 ds ≤ i(0) 2 e -µit + C(1 -e -µit ).

Continuity of the operator R

Let (v 1 , i 1 ) → (W 1 , I 1 , R 1 ) → (v 1 , i 1 ) and (v 2 , i 2 ) → (W 2 , I 2 , R 2 ) → (v 2 , i 2
) be two solutions of the systems ( 7)-( 8), with the same initial conditions (v 0 = v 0 , i 0 = i 0 , W 0 , I 0 , R 0 , v 0 , i 0 ).

Define v = v 2 -v 1 , î = i 2 -i 1 , Ŵ = W 2 -W 1 , Î = I 2 -I 1 , R = R 2 -R 1 , v = v 2 -v 1 , î = i 2 -i 1 .
Lemma 3 There exists a positive continuous function C(t) such that

t 0 v 2 + î 2 ds ≤ C(t) t 0 v 2 + î 2 ds . (9)
Proof: The differences v, î, Ŵ , Î, R solve the system

∂ t v -d v ∆ î + µ v v + α 4 W 2 v + α 4 v1 Ŵ -α v Î = 0 , ∂ t î -d i ∆ î + µ i î + α 3 W 2 î + α 3 i 1 Ŵ -α i Î = 0 , ∂ t Ŵ -d∆ Ŵ + i 2 Ŵ + W 1 î + v 2 Ŵ + W 1 v = 0 . ∂ t Î -d∆ Î + µ I Î + W 2 v + v 1 Ŵ = 0 , ∂ t R -d∆ R + i 2 Ŵ + W 1 î = 0 .
with homogeneous initial and boundary conditions.

Multiply the equation for Î by Î and integrate by parts to obtain 1 2

d dt Î 2 + d ∇ Î 2 + µ I Î 2 = - Ω v W 2 Î dΩ + Ω v 1 Ŵ Î dΩ .
Noticing that all the functions are in L ∞ , we deduce that there exists C > 0, depending on the initial conditions and the model coefficients, such that

d dt Î 2 + 2d ∇ Î 2 ≤ C v 2 + Ŵ 2 + Î 2 . ( 10
)
Using the same argument we can show that

d dt Ŵ 2 + 2d ∇ Ŵ 2 ≤ C î 2 + v 2 + Ŵ 2 , (11)
d dt R 2 + 2d ∇ R 2 ≤ C î 2 + Ŵ 2 .
From (10)-(11) and d ≥ 0, it follows that

d dt Ŵ 2 + Î 2 + R 2 ≤ C î 2 + v 2 + Ŵ 2 + Î 2 .
Integrating over [0, t] we obtain that

Ŵ 2 + Î 2 + R 2 ≤ C t 0 î 2 + v 2 ds + C t 0 Ŵ 2 + Î 2 ds .
Applying Gronwall's Lemma it follows that there exists a positive continuous function C(t) such that

Ŵ 2 + Î 2 + R 2 ≤ C(t) t 0 î 2 + v 2 ds . (12)
On the other hand, multiply the equation for î by î and integrate by parts to obtain 1 2

d dt î 2 + d i ∇ î 2 + µ i î 2 + α 3 Ω W 2 î2 dΩ = -α 3 Ω i 1 Ŵî dΩ + α i Ω Îî dΩ .
Proceeding as before we can show that there exists a constant C 1 > 0 such that

d dt î 2 + 2d i ∇ î 2 ≤ C 1 î 2 + Ŵ 2 + Î 2 .
Analogously, we can show that

d dt v 2 + 2d v ∇v 2 ≤ C 1 v 2 + Ŵ 2 + Î 2 .
Using (12) we obtain that

d dt î 2 + v 2 ≤ C 1 î 2 + v 2 + C(t) t 0 î 2 + v 2 ds .
Integrating over [0, t] yields

î 2 + v 2 ≤ C 1 t 0 î 2 + v 2 ds + C(t) t 0 î 2 + v 2 ds .
Applying Gronwall's Lemma, it follows that

î 2 + ŵ 2 ≤ C(t) t 0 î 2 + v 2 ds .
Integrating again over [0, t] we obtain [START_REF] Taylor | Partial differential equations III[END_REF].

5 Proof of the theorems

Proof of Theorem 1

Lemma 4 Fix a positive time T > 0. Then:

1. K is a convex closed subset of L 2 0, T ; L 2 (Ω) 2 . 2. R : K → K. 3. R : L 2 0, T ; L 2 (Ω) 2 → L 2 0, T ; L 2 (Ω) 2 is continuous. 4. R(K) is relatively compact in L 2 0, T ; L 2 (Ω) 2 .
Proof:

1. By construction, K is convex and closed.

2. If i ≥ 0 and v ≥ 0 then from Lemma 1 it follows that i ≥ 0 and v ≥ 0.

3. By Lemma 3 the operator R is continuous.

4. We will use Aubin's compactness theorem (see Theorem 5.1 in Lions [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites nonlinéaires[END_REF], Section 5.5, pp. 57-64, and Tartar [START_REF] Tartar | An introduction to Navier-Stokes equation and Oceanography[END_REF], Chapter 24, pp. 137-141). Suppose that the sequence

{(v n , i n )} is uniformly bounded in L 2 0, T ; L 2 (Ω) 2 .
Then, by the continuity of R, the sequence {R[(v n , i n )] = (v n , i n )} is also uniformly bounded in L 2 0, T ; L 2 (Ω) , and the estimates in Lemma 2 imply that

{(v n , i n )} is uniformly bounded in L 2 0, T ; H 1 (Ω) 2 . Furthermore, the sequence of derivatives {(∂ t v n , ∂ t i n )} is uniformly bounded in L 2 0, T ; H -1 (Ω)
2 , because the expressions

∂ t v n = ∇ • (d v ∇v n ) -µ v v n + α v I n -α 4 v n W n , ∂ t i n = ∇ • (d i ∇i n ) -µ i i n + α i I n -α 3 i n W n ,
define two uniformly bounded sequences of distributions in H 1 (Ω). Therefore, applying Aubin's theorem to the spaces

H 1 (Ω) 2 ⊂ L 2 (Ω) 2 ⊂ H -1 (Ω) 2 , we obtain that the sequence {(v n , i n )} is relatively compact in L 2 0, T ; L 2 (Ω) 2 .
We can now conclude the existence of solutions of the reaction-diffusion system ( 2)-( 3):

From Lemma 4 we see that the operator R satisfies Schauder's Fixed Point Theorem (Corollary B.3, p. 262 in Taylor [?]). Therefore, the system (2)-( 3) has solutions W, I, R, v, i. Moreover, from the fact that v and i coincide with v and i , respectively, we can derive two consequences. First, from Lemma 1 it follows that W, I, R, v, i are non-negative and bounded. Second, from Lemma 2 we have that W, I, R, v, i belong to L 2 0, T ; H 1 (Ω) 2 .

Lemma 5 The solutions W (x, t), I(x, t), R(x, t), v(x, t) and i(x, t) of the reaction-diffusion system (2)-( 3) are unique and depend continuously on the initial data (4).

Proof: Notice that v and i coincide with v and i , respectively. Bearing this in mind, and repeating the arguments in Lemma 3, we can show that there is a positive constant C such that

d dt Ŵ 2 + Î 2 + R 2 + v 2 + î 2 ≤ C Ŵ 2 + Î 2 + R 2 + v 2 + î 2 .
Therefore, Gronwall's lemma implies that there exists a positive continuous function C(t) such that

Ŵ 2 + Î 2 + R 2 + v 2 + î 2 ≤ C(t) Ŵ0 2 + Î0 2 + R0 2 + v0 2 + î0 2 .

Proof of Theorem 2

So far we have only considered the reaction-diffusion system (2)-(3). However, for the hybrid system ( 5)-( 6) the same results hold. Indeed, on one hand, in Lemmas 2-5 we have only used that d ≥ 0. On the other hand, we can use the integral representations of W , I and R in order to deduce that they are bounded and non-negative, which proves Lemma 1 in the hybrid case.

Proof of Theorem 3

For any 0

< d ≤ d 0 let (W d , I d , R d , v d , i d
) be a weak solution of ( 2)-( 3) and let d → 0.

First, from Lemma 2 it follows that the sequence W d is bounded in

L 2 0, T ; H 1 (Ω) ∩ L ∞ (Ω × (0, T )) ,
which implies that a subsequence, still denoted W d , converges weakly in L 2 0, T ; H 1 (Ω) and weaklyin L ∞ (Ω × (0, T )) to W 0 . This result is also valid for the sequences

I d , R d v d and i d .
Second, using classical estimates of type

T 0 Ω |i d W d -i 0 W 0 | dΩdt ≤ C 1 i d -i 0 + C 2 W d -W 0
it can be shown that the sequence ∂ t W d is uniformly bounded in L 2 0, T ; H -1 (Ω) . In consequence, applying Aubin's compactness theorem yields the strong convergence i d → i 0 in L 2 0, T ; L 2 (Ω) . This result is also valid for the sequences I d , R d v d and i 0 .

Finally, we obtain that a subsequence (W d , I d , R d , v d , i d ) of weak solutions of ( 2)-( 3) converges strongly in L 2 0, T ; L 2 (Ω) , weakly in L 2 0, T ; H 1 (Ω) and weakly-in L ∞ (Ω × (0, T )). Since the limit (W 0 , I 0 , R 0 , v 0 , i 0 ) is by construction a weak solution of ( 5)-( 6), the uniqueness of the hybrid system ( 5)-(6) implies that the limit is the same for any converging subsequence.

Proof of Theorem 4

Lemma 6 The solutions W, I, R, v, i of the RD system (2)-(3) are globally-defined and belong to L ∞ (Ω × (0, ∞)).

Proof:

1. Define N (x, t) := W (x, t) + I(x, t) + R(x, t). Then N satisfies ∂ t N -d∆N ≤ 0 in Ω × [0, T ], ∇N • n = 0 on Γ × [0, T ].
Therefore, Maximum Principle provides

N (x, t) ≤ N 0 ∞ ∀(x, t) ∈ Ω × [0, T ].
Moreover, since this bound is independent of t and T , it follows that N (x, t) is defined for all t ∈ R.

In consequence, the positivity of W, I, R implies that these three functions exist for all times and that they are uniformly bounded, i.e.,

W (x, t), I(x, t), R(x, t) ≤ N 0 ∞ ∀(x, t) ∈ Ω × [0, ∞).

Proof of Theorem 5

Lemma 8 If W, I, R, i, v are non-negative, steady-state solutions of the hybrid system ( 5)-( 6) then

I(x) ≡ 0, v(x) ≡ 0, i(x) ≡ 0.
Proof: The steady-state-solutions of system ( 5)-( 6) solve

           0 = -iW -vW , 0 = -µ I I + vW , 0 = iW , -d v ∆v = -µ v v + α v I -α 4 vW , -d i ∆i = -µ i i + α i I -α 3 iW , (15) with boundary conditions ∇v • n(σ) = 0 on Γ, ∇i • n(σ) = 0 on Γ. ( 16
)
From the first tree equations in (15)-( 16) it follows immediately that

iW = vW = µ I I = 0.
Since µ i > 0 then I = 0. Moreover, v and i solve

-d v ∆v = -µ v v, -d i ∆i = -µ i i.
In consequence, since µ v > 0 and µ i > 0 it follows that v = i = 0. It is worth to mention that in this case we have no restriction on W and R.

Lemma 9 Suppose that the initial conditions belong to L ∞ (Ω). For any x ∈ Ω, the solutions of the hybrid system ( 5)-( 6) are integrable in the following sense:

1.

lim t→∞ t 0 I(x, s) ds < ∞ and lim t→∞ t 0 Ω I(x, s) dΩ ds < ∞.
2.

lim t→∞ t 0 v(x, s)W (x, s) ds < ∞ and lim t→∞ t 0 Ω v(x, s)W (x, s) dΩ ds < ∞.
3.

lim t→∞ t 0 i(x, s)W (x, s) ds < ∞ and lim t→∞ t 0 Ω i(x, s)W (x, s) dΩ ds < ∞.
Proof:

1. Adding the equations of W, I, R we obtain that

∂ t W + ∂ t I + ∂ t R + µ I I ≤ 0.
Integrating over [0, t] it follows that, for any (x, t) ∈ Ω × [0, ∞) we have

W (x, t) + I(x, t) + R(x, t) + µ I t 0 I(x, s)ds ≤ W (x, 0) + I(x, 0) + R(x, 0). ( 17
)
Therefore, from (17) and using uniform boundedness of the initial conditions we can deduce that lim Lemma 10 The solutions of the hybrid system satisfy the following properties:

< t 0 < t 1 < • • • < t n → ∞
and define f n (x) := I(x, t n ).

We have already proven that For the last statement, if we integrate the equation for I over Ω × [0, t] we obtain

I(t) 2 + 2µ I t 0 Ω I 2 (x, s) dΩ ds = I(0) 2 + 2 I L ∞ (Ω×(0,t)) t 0 Ω v(x, s)W (x, s) dΩ ds.
Since the right-hand side converges as t → ∞ then the left converges as well, which implies that

lim t→∞ t 0 Ω I 2 (x, s) dΩ ds < ∞.
2. Multiplying the equation for v by v and integrating by parts yields 1 2

d dt v(t) 2 + d v ∇v(t) 2 + µ v v(t) 2 = α v Ω Iv dΩ -α 4 Ω v 2 W dΩ.
Integrating over [0, t] it can be shown that there is a constant C > 0, independent of t, such that

v(t) 2 + 2d v t 0 ∇v(s) 2 ds + µ v t 0 v(s) 2 ds ≤ v(0) 2 + C t 0 Ω I 2 (x, s) dΩ ds.
But from the previous results it follows that the right-hand side converges. In consequence, the left-hand side converges, i.e., lim

t→∞ t 0 v(s) 2 + ∇v(s) 2 ds < ∞.
The proof for i is exactly the same as for v.

3. Let ϕ(x, t) be a test function in L 1 0, ∞; H 1 (Ω) . Calculating the dual product in L 1 0, t; H 1 (Ω) yields

∂ t v, ϕ = -d v t 0 Ω ∇v • ∇φ dΩ dt -µ v t 0 Ω
vϕ dΩ dt

+α v t 0 Ω Iϕ dΩ dt -α 4 t 0 Ω
vW ϕ dΩ dt .

Using the uniform boundedness and integrability of I and W it can be shown that there exist three positive constants C 1 , C 2 and C 3 , independent of t, such that

| ∂ t v, ϕ | ≤ C 1 + C 2 ϕ L 1 (0,t;H 1 (Ω)) + C 3 v L 1 (0,t;H 1 (Ω)) .
This implies that | ∂ t v, ϕ | is bounded for any test function in L 1 0, t; H 1 (Ω) with norm less than one, uniformly in t. Therefore,

∂ t v is a bounded distribution in L 1 0, ∞; H 1 (Ω) .
The proof for i is exactly the same as for v.

Lemma 11 Suppose that the initial conditions belong to L ∞ (Ω). For any x ∈ Ω fixed, the solutions of the hybrid system ( 5)-( 6) satisfy the following properties: From the previous results and the uniform boundedness of W , I and v it follows that

lim t→∞ t 0 Ω v(x, t) dΩ dt < ∞. Since ∂ t v is uniformly bounded it follows that lim t→∞ Ω v(x, t) dΩ = 0. ( 18
)
Now let ε > 0 and choose x 0 ∈ Ω. Define

B(ε, x 0 ) = {x ∈ Ω : |x -x 0 | < ε}.
From (18) there exists T (ε) > 0 such that Ω v(x, t) dΩ < ε|B(ε, x 0 )| for all t ≥ T (ε).

In consequence,

Ω v(x, t) dΩ ≤ Ω v(x, t) dΩ < ε|B(ε, x 0 )|, which implies that 1 |B(ε, x 0 )| Ω v(x, t) dΩ ≤ ε.
Therefore, since v(x, t) is continuous, we have

v(x 0 , t) = lim |B(ε,x0)|→0 1 |B(ε, x 0 )| Ω v(x, t) dΩ ≤ ε, which implies that lim t→∞ v(x 0 , t) ≤ ε.
Since ε > 0 and x 0 ∈ Ω are arbitrary it follows that lim t→∞ v(x, t) = 0 for all x ∈ Ω.

The proof for i is exactly the same.

Define

C ∞ := lim t→∞ t 0 Ω v(x, s) dΩ ds.
From the previous results we have that C ∞ is well-defined, finite and positive. Now, for any n ∈ N define

A n := x ∈ Ω : lim t→∞ t 0 v(x, s) ds > nC ∞ .
On the one hand,

lim t→∞ t 0 An v(x, s) dΩ ds > nC ∞ |A n | .
On the other hand,

lim t→∞ t 0 An v(x, s) dΩ ds ≤ lim t→∞ t 0 Ω v(x, s) dΩ ds = C ∞ .
Therefore,

|A n | < 1 n . (19)
Moreover, since In consequence, the result follows immediately from Lemma 11.

A 1 ⊃ A 2 ⊃ • • • ⊃ A n ⊃ • • • it follows that

Proof of Theorem 6

Lemma 13 If v ∞ (x) is a steady-state solution of the hybrid system (5)-( 6) and µ v = 0 then ∇v ∞ = 0.

Proof: Following the proof of Lemma 8 it follows that

d v ∆v ∞ = 0.
Therefore, integrating by parts yields ∇v ∞ = 0.

Lemma 14 If v is solution of the hybrid system (5)-( 6) and µ v = 0 then

Discussion

We proved that both the reaction-diffusion and the hybrid models are well-posed problems, i.e., they have global unique solutions, which are non-negative, bounded and depend continuously on the initial data. We also showed that, when d → 0, the solution of the reaction-diffusion model coverges to the solution of the hybrid model.

We provided several asymptotic estimates for the solutions of the hybrid model. First, we proved tat the solutions are uniformly bounded and integrable over Ω × (0, ∞). Second, we showed that virions and wild-type cells cannot coexist, because the product of their concentrations tends to zero as t → ∞. Third, we proved that if µ v > 0 then the virus concentration tends to zero as t → ∞ and W tends to a non-zero, non-homogeneous limit. Fourth, if µ v = 0 and α v ≥ α 4 µ I then the global virus concentration is bigger than the original concentration.

One striking result for the hybrid model is the global stability of the steady-state solutions. Indeed, we characterized the steady-state solutions and showed that they coincide with the limits of the corresponding time-dependent solutions. Indeed, for µ v > 0, lim Finally, in the numerical simulations we found that the spacial structure (i.e. the diffusion) plays a crucial role in the proliferation of virions. Indeed, d v has a positive effect when µ v = 0 and a negative effect when µ v > 0. Therefore, the spacial structure permits the existence of richer patterns than in the original ODE system of P. Getto et al [START_REF] Getto | Modelling and analysis of dynamics of viral infection of cells and of interferon resistance[END_REF], which confirms what we have conjectured at the beginning of the project.

Figure 1 :

 1 Figure 1: dv = 0.001, µv = 0.2 and αv = 0.8.

Figure 2 :

 2 Figure 2: dv = 0.01, µv = 0.2 and αv = 0.8.

Figure 3 :

 3 Figure 3: dv = 0.1, µv = 0.2 and αv = 0.8.

Figure 4 :

 4 Figure 4: dv = 1, µv = 0 and αv = 2.

Figure 5 :

 5 Figure 5: dv = 0.001, µv = 0 and αv = 2.

Figure 6 :

 6 Figure 6: dv = 0.0001, µv = 0 and αv = 2.

3 .

 3 s) dΩ ds < ∞. 2. Integrating over [0, t] the equation for I yields I(x, t) -I(x, 0) + t 0 I(x, s) ds = t 0 v(x, s)W (x, s) ds. Therefore, using the previous result on the integrability of I we obtain that lim t→∞ t 0 v(x, s)W (x, s) ds < ∞ and lim t→∞ t 0 Ω v(x, s)W (x, s) dΩ ds < ∞. Integrating over [0, t] the equation for W yields W (x, t) -W (x, 0) + t 0 vW (x, s) ds = t 0 i(x, s)W (x, s) ds. Since the left-hand side has a limit so does the right-hand side, i.e., lim t→∞ t 0 i(x, s)W (x, s) ds < ∞ and lim t→∞ t 0 Ω i(x, s)W (x, s) dΩ ds < ∞.

  lim n→∞ f n (x) = 0 and 0 ≤ f n (x) ≤ N 0 (x), where N 0 (x) := W (x, 0) + I(x, 0) + R(x, 0) ∈ L ∞ (Ω).Therefore, we can apply Lebesgue's Dominated Convergence Theorem to obtain lim for any arbitrary sequence (t n ) n∈N then lim t→∞ Ω I 2 (x, t) dΩ = 0.

  , s)W (x, s) dΩ dt.

  In consequence, from (19) we obtain that |A ∞ | = 0, which implies that lim t→∞ t 0 v(x, s) ds < ∞ a.e. in Ω.Finally, since t → v(x, t) is continuous we conclude thatlim t→∞ t 0 v(x, s) ds < ∞ for all x ∈ Ω.The proof for i is the same.Lemma 12 For any x ∈ Ω, W 0 (x) > 0 if and only if lim t→∞ W (x, t) > 0.Proof: Integrating over [0, t] the equation for W we obtain W (x, t) = W 0 (x) × exp -

  x) dΩ + (α v -α 4 µ I ) Integrating over Ω × [0, t] the equation for v and using Fubini's Theorem yields Ω s)W (x, s) dΩ dt. Now recall that v satisfies homogeneous Neumann boundary conditions. Therefore, t 0 Ω ∆v(x, t) dΩ dt = t 0 ∂Ω ∇v(x, t) • n dS dt = 0,

  t→∞ I(x, t) = 0 = I ∞ (x) , lim t→∞ v(x, t) = 0 = v ∞ (x) , lim t→∞ i(x, t) = 0 = i ∞ (x) ,whilst W (x, t) has an explicit limit, W ∞ (x) = W 0 (x) exp -

  1. ∂ t I is uniformly bounded for t ∈ [0, ∞). Moreover, for any x ∈ Ω, ∂ t v and ∂ t i are bounded distributions in L 1 0, ∞; H 1 (Ω) .

	lim t→∞ lim t→∞ 1. Since the functions I, v and W are uniformly bounded in Ω × (0, ∞), from the equation for I it I(x, t) = 0, lim t→∞ Ω I 2 (x, t) dΩ < ∞ and lim t→∞ t 0 Ω I 2 (x, s) dΩ ds < ∞. 2. t 0 v(s) 2 + ∇v(s) 2 ds < ∞ and lim t→∞ t 0 i(s) 2 + ∇i(s) 2 ds < ∞. follows that ∂ t I is uniformly bounded as well. Therefore, using lim t→∞ t 0 I(x, s) ds < ∞ we can infer that lim t→∞ I(x, t) = 0 for all x ∈ Ω. Now consider a sequence 3. Proof: 0

  On one hand, using the previous results it follows that the right-hand side is uniformly bounded in t. In consequence, On the other hand,∂ t (vW ) = v∂ t W + W ∂ t vis uniformly bounded because v, W , ∂ t v and ∂ t W are all uniformly bounded. In consequence,

				t
			lim t→∞	0	v(x, t)W (x, t)ds < ∞.
			lim t→∞	v(x, t)W (x, t) = 0.
	2. Integrating over Ω × [0, t] the equation for v and using Fubini's Theorem yields
				t		t
					∆v(x, t) dΩ dt -µ v	v(x, t) dΩ dt
				0	Ω	0	Ω
					t	t
			+α v	I(x, s) dΩ dt + α 4
					0	Ω	0	Ω
	1.
	lim t→∞	v(x, t)W (x, t) = 0 and lim t→∞	i(x, t)W (x, t) = 0,
	2.
			lim t→∞	v(x, t) = lim t→∞	i(x, t) = 0.
	3.	t				t
	lim t→∞	0	v(x, s) ds < ∞ and lim t→∞	0	i(x, s) ds < ∞ .
	Proof:
	1. Integrating the equation for I over [0, t] we obtain

t 0 v(x, s)W (x, s)dt = I(x, t) -I(x, 0) + µ I t 0 I(x, s)ds .

Ω

v(x, t) dΩ -Ω v(x, 0) dΩ = v(x, s)W (x,

s) dΩ dt. Now recall that v satisfies homogeneous Neumann boundary conditions. Therefore, t 0 Ω ∆v(x, t) dΩ dt = t 0 ∂Ω ∇v(x, t) • n dS dt = 0, which implies that

2.

From Lemma 1 we obtain that 0 ≤ v(x, t) ≤ γ 4 and 0 ≤ i(x, t) ≤ γ 5 for all (x, t) ∈ Ω × (0, T ), where

Again, since the bounds are independent of t and T , the solutions v and i exist for all times.

Lemma 7 If W, I, R, i, v are non-negative, steady-state solutions of the RD system (2)-(3) then

Proof: If W, I, R, i, v are non-negative, steady-state solutions of system (2)-( 3), then

with boundary conditions

Multiplying the equation for W by W and integrating by parts we obtain

Therefore iW = vW = 0 and W (x) = W 0 constant. Using same argument for I yields

which implies that I(x) ≡ 0. For R we obtain

and therefore, R(x) = R 0 is constant, whilst for v, i we obtain

and, in consequence, v(x) = i(x) ≡ 0. . which implies that

On the other hand, integrating the equation for I yields

In consequence,

In particular, if v 0 ≡ 0 then v ∞ ≡ 0.

Proof: If α v ≥ α 4 µ I then from Lemma 14 we have Recall that v(x, t) is uniformly bounded. Therefore, applying Fatou's lemma we obtain