G Bocharov

A Meyerhans

N Bessonov

S Trofimchuk

V Volpert
email: volpert@math.univ-lyon1.fr

Modelling the dynamics of virus infection and immune response in space and time

Keywords: infection spreading, tissue, immune response, mathematical model, reaction-diffusion equations, time-delay

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Acute and chronic virus infections remain to be a problem of major public concern worldwide. Defence of the host organism against an infectious challenge is a primary function of the immune system associated with generation of potent responses [START_REF] Qi | Spatiotemporal basis of innate and adaptive immunity in secondary lymphoid tissue[END_REF]. Understanding the processes and parameters that underline fate decision between resolving and chronic dynamics of the virus infections such as CMV, HBV, HCV, HIV is a challenge. One of the major goals of systems approaches to immunology is to quantitatively measure complex multiscale spatial and temporal interactions, permitting development of computational models that can be used to predict the responses to perturbations [START_REF] Subramanian | Network representations of immune system complexity[END_REF]. In addition, they feed conceptual developments in theory of the immune system which due to its capability to perform many complex "computations" in a highly parallel and distribute fashion, is considered as a source of inspiration for developing intelligent methodologies for complex problem solving [START_REF]Artificial Immune Systems and Their Applications[END_REF].

Although the mathematical models, like experimental animal systems [START_REF] Vrisekoop | Life and death as a T lymphocyte: from immune protection to HIV pathogenesis[END_REF], can not entirely reproduce the complex etiology of virus infections such as HIV in hu- January 9, 2017 20:38 The International Journal of Parallel, Emergent and Distributed Systems imm-sys3 2 Bocharov, Meyerhans, Bessonov, Trofimchuk, Volpert mans or experimental LCMV infection, the modelling results might inform studies more difficult to undertake in humans. In particular, the computational models of antiviral immune response can provide insights into the regulation of the dynamics processes underlying various patterns in disease progression. Except for few studies, mathematical models of the immune response against virus infections conventionally consider the infected whole organism as a single homogenous compartment [START_REF] Cappuccio | Multiscale modelling in immunology: a review[END_REF][START_REF] Bocharov | Spatiotemporal dynamics of virus infection spreading in tissues[END_REF]. Thus they do not take into account that the dynamics of infection spread differs between tissues, organs and blood. A more detailed description however should include this spatial heterogeneity with respect to virus propagation and immune response development [START_REF] Bouchnita | Hybrid approach to model the spatial regulation of T cell responses[END_REF], [START_REF] Bouchnita | Towards a Multiscale Model of Acute HIV Infection[END_REF].

Development of a viral infection (e.g., HIV, SIV, LCMV) in tissues such as lymph nodes or spleen occurs due to virus spreading that may be either direct between neighboring cells (cell-to cell model) or by means of a random motion of the free virus in the extracellular matrix. The infected cells are the sources of virus production. Reaction-diffusion equations represent an appropriate model for the description of virus transport and production, and of the spatiotemporal dynamics immune response that eliminates infected cells.

Interaction between the virus and the host organism includes many cellular and humoral components of the innate and adaptive arms of the immune system [START_REF] Zinkernagel | On immunity against infections and vaccines: credo 2004[END_REF]. Virus-specific cytotoxic T lymphocytes (CTL) play a key role in the dynamics and outcome of virus infections [START_REF] Moskophidis | Virus persistence in acutely infected immunocompetent mice by exhaustion of antiviral cytotoxic effector T cells[END_REF][START_REF] Mcmichael | The immune response during acute HIV-1 infection: clues for vaccine development[END_REF] . Importantly, that CTL temporal dynamics and viral growth dynamics can under many circumstances be successfully captured using simple conceptual views of the regulation based on predator-prey systems [START_REF] Nowak | Virus Dynamics[END_REF], there are other phenomena, such as CTL exhaustion, which do not fit in this simple predator-prey interpretation and demand more complex conceptual models (e.g., the balance of growth and differentiation models) to capture reality [START_REF] Grossman | Dynamic tuning of lymphocytes: physiological basis, mechanisms, and function[END_REF]. Since long ago, it was observed that virus infections have a strong impact on the immune system [START_REF] Notkins | Effect of Virus Infections on the Function of the Immune System[END_REF].

The experimental LCMV infection in mice provided a solid evidence that the intensity of the CTL response depends on the viral load and growth kinetics in a nonlinear manner [START_REF] Moskophidis | Virus persistence in acutely infected immunocompetent mice by exhaustion of antiviral cytotoxic effector T cells[END_REF][START_REF] Bocharov | Understanding Experimental LCMV Infection of Mice: The Role of Mathematical Models[END_REF]. A small infection level will stimulate an immune response while a sufficiently high infection level will downregulate an immune response by immune cell death and/or decrease of immune cell proliferation. Moreover the action of the immune response is delayed in time since, after the stimulation by infection, immune cells need some time to proliferate and differentiate. Thus, infection spreading in the tissue can be modelled by reaction-diffusion equations with delay.

In our previous work [START_REF] Bocharov | Spatiotemporal dynamics of virus infection spreading in tissues[END_REF] we studied the spread of a viral infection in tissue with a single reaction-diffusion equation for the virus concentration. The concentration of the immune cells was considered as a function of the infection level. With this model different regimes of infection spreading were identified. Depending on the intensity of the immune response and the initial viral load, we observed a) infection propagation with a low viral load established behind the infection front, b) a high infection level behind the infection front, and c) two consecutive infection fronts, the first one with a low infection level and the second one behind it with a high infection level.

In this work here we now consider a more detailed model of infection spreading in tissue with a reaction-diffusion system of two equations, one for the concentration of virus and one for that of immune cells (CTLs). The model is formulated in the next section. We discuss the corresponding ordinary differential system of equations in Section 3. The results of numerical simulations of the reaction-diffusion system with delay are presented in Section 4. Finally, we discuss the results and implications of January 9, 2017 20:38 The International Journal of Parallel, Emergent and Distributed Systems imm-sys3 Parallel, Emergent and Distributed Systems 3 our analysis for understanding the pathogenesis of virus infections in Section 5.

Model

The biological scheme of the model is presented in Figure 1. In the model we consider multiple non-linear feedback regulations of the population dynamics of virus infection and the antiviral immune responses. According to this general scheme, the viruses and immune cells compete for their survival. The outcome of the competition results from the race between the virus infection spread on one side and the clonal CTL expansion followed by migration to the site of infection to eliminate the virus sources. Depending on the viral load and growth kinetics, the efficacy of the CTL response in determined by a balance of proliferation, anergy induction and activation-induced cell death by apoptosis. The later two take into account the tuning of lymphocytes [START_REF] Grossman | Dynamic tuning of lymphocytes: physiological basis, mechanisms, and function[END_REF], a key feature of the immune system.

Two PDEs model

We consider the following system of equations for the concentration of virus v and of immune cells c in lymphoid organs such as the spleen or lymph nodes, viewed as 1D structures:

∂v ∂t = D 1 ∂ 2 v ∂x 2 + kv(1 -v) -σvc, ( 1
)
∂c ∂t = D 2 ∂ 2 c ∂x 2 + (c 0 + ϕ(v τ )c)(1 -c) -ψ(v µ )c. (2)
Here v and c are the normalized dimensionless concentrations of virus and immune cells, respectively, ψ(v µ ) will be specified below. The first term in the right-hand side of equation ( 1) describes the virus diffusion, the second term its production and the last term its elimination by the immune cells. The parameter D 1 is the diffusion coefficient of virus and k 1 stands for the replication rate constant. Similarly, the first term in the right-hand side of the second equation describes random motion of immune cells, the second term gives the rate of the appearance of immune cells migrating to the peripheral organs from the bone marrow and thymus and proliferating there and the last term parameterizes the rate of their death. The second term includes cell influx due to differentiation of less mature precursor cells and cell birth due to antiviral immune response. The latter, depends on the viral load in a non-linear bell-shaped manner. The last term for cell death includes natural (constant) cell death and their apoptosis induced by infection. Both of them are included in the function ψ(v µ ). Parameters D 1 , D 2 , k and σ are positive numbers, the functions ϕ(v) and ψ(v) are assumed to be positive and sufficiently smooth. The delay values in ϕ(v) and ψ(v) represent the net duration of cell proliferation and differentiation to effector cells and the time needed for programming of activated CTLs to apoptosis, accordingly [START_REF] Bocharov | Modelling the dynamics of LCMV infection in mice: conventional and exhaustive CTL responses[END_REF].

v τ = v(x, t -τ ), v µ = v(x, t -µ), the functions ϕ(v τ ) and
In what follows we will consider the function

f (v) = 1 -ψ(v)/ϕ(v).
We will suppose that f (0) > 0, it is increasing for sufficiently small v and decreasing for v large enough. Such properties of this function corresponds to the biological observation that small viral load promotes proliferation of immune cells while large level of infection decreases cell proliferation and increases their death. This function can remain positive for all v or it can become negative for v sufficiently large. A typical form of this function is shown in Figure 2.

We note that the number of immune cells increases enormously i.e. by 10 3 -fold, via clonal expansion in the course of immune response to infection. Therefore, we will neglect the cell influx term in comparison with the effect of clonal proliferation on the population size and we will set c 0 = 0. This assumption allows us to simplify the model analysis.

Reduction to a single PDE model

If we consider the functions ϕ and ψ in the form

ϕ(v) = Rϕ 0 (v), ψ(v) = Rψ 0 (v),
where R is a large parameter, then we can divide the second equation by R and formally pass to the limit R → +∞. Then we get

c = f (v), where f (v) = 1 - ψ(v) ϕ(v) .
If τ = µ, then substituting this expression into equation (1), we get the equation

∂v ∂t = D 1 ∂ 2 v ∂x 2 + kv(1 -v) -σf (v τ )v (3)
introduced and studied in [START_REF] Bocharov | Spatiotemporal dynamics of virus infection spreading in tissues[END_REF]. As it is mentioned in the introduction, there are three main regimes of the infection spreading in this case: with low level infection level behind the infection front, with high level infection level, and with low level followed by high level infection. The choice between them is determined by the intensity of immune response (function f (v)) and by the initial infection distribution in the tissue.

Existence and stability of stationary points

Two ODEs version of the model

We begin to study system (1), ( 2) with the stationary points of the ODE system

dv dt = kv(1 -v) -vc, ( 4
)
dc dt = ϕ(v)c(1 -c) -ψ(v)c, ( 5
)
where we set σ = 1. The stationary points of this system are determined by the equations Some typical examples of nullclines, stationary points and trajectories of this system are shown in Figure 2. There are 6 stationary points in the left figure, P 1 = (0, 0), P 2 = (0, f (0)),

v = 0 or c = k(1 -v); c = 0 or c = f (v). (6)
P i = (v i , c i ), i = 3, 4, 5, 6, where 0 < v 3 < v 4 < v 5 < v 6 = 1, 0 = c 6 < c 5 < c 4 < c 3 .
Stability of these points are determined by the eigenvalues of the matrix

A = ( k(1 -2v) -c -v ϕ ′ (v)c(1 -c) -ψ ′ (v)c ϕ(v)(1 -2c) -ψ(v)
)

. It can be easily verified that P 1 and P 6 are saddle points, P 2 is also a saddle if k > f (0) (Figure 2, left) and it is a stable node if k < f (0) (Figure 2, right). The other points with positive coordinates are stabile if f ′ (v i ) > -k and unstable if the inequality is opposite.

Depending on the values of parameters, the trajectory starting in the vicinity of the unstable point P 1 can have one of the following behaviors.

1. If f (0) < k, then the point P 3 (0, f (0)) is stable (Figure 2, left). Since v 3 > 0, then virus is not eliminated. This case corresponds to chronic infection with a low viral load. The level of virus infection can monotonically increase in time converging to its limiting value v 3 or it can pass over its maximal value and decrease towards v 3 .

2. If f (0) > k, then the point P 2 (0, f (0)) is stable, and the trajectory can converge to it for large time (Figure 2, right). This case corresponds to complete virus elimination after acute infection. The memory immune cells remain in the body providing rapid immune response in the case of further infections.

3. Another stable stationary point is P 5 = (v 5 , c 5 ), where f ′ (v) > -1. The corresponding value of v is large and of c low. This is the case of infection induced immunodeficiency characterized by high viral load and partial exhaustion of CTLs. 2, right), then the stable stationary point v = 1, c = 0 corresponds to complete exhaustion of immune cells.

4. If f (v) < 0 for v = 1 (Figure

Two Delay Differential Equations version of the model

In this section we study stability of the same stationary points of the ODE system with the delays introduced into the corresponding terms as in the original PDE model ( 1)-(2):

dv dt = kv(1 -v) -vc, ( 7
)
dc dt = ϕ(v τ )c(1 -c) -ψ(v µ )c. (8)
Let us note that a similar system with particular form of the functions ϕ(v) and ψ(v) was studied in [START_REF] Bocharov | Modelling the dynamics of LCMV infection in mice: conventional and exhaustive CTL responses[END_REF]. In order to analyze below behavior of solutions of the reaction-diffusion system, we begin with the study stability of the stationary solution (v * , c * ) with respect to this DDEs system with two delays. After linearization about this stationary solution, we obtain the equation with respect to the eigenvalue λ:

det(A(λ) -λE) = 0, (9)
where

A(λ) = ( -kv * , -v * c * g(λ) , -c * ϕ(v * ) ) , g(λ) = ϕ ′ (v * )(1 -c * )e -λτ -ψ ′ (v * )e -λµ .
We recall that

c * = k(1 -v * ), ϕ(v * )(1 -c * ) -ψ(v * ) = 0. ( 10
)
Equation ( 9) can be written as follows:

ϕ(v) (kv + λ)(k(1 -v)ϕ(v) + λ) kv(1 -v) = ϕ(v)ψ ′ (v)e -λµ -ϕ ′ (v)ψ(v)e -λτ (11)
(v * is replaced here by v for simplicity of notation). In order to find the onset of oscillatory instability, we set λ = iξ. Equating the real and imaginary parts, we get: 12)

ϕ(v) kv(1 -v) (k 2 v(1 -v)ϕ(v) -ξ 2 ) = ϕ(v)ψ ′ (v) cos(ξµ) -ϕ ′ (v)ψ(v) cos(ξτ ), (
ϕ(v) v(1 -v) (v + (1 -v)ϕ(v))ξ = -ϕ(v)ψ ′ (v) sin(ξµ) + ϕ ′ (v)ψ(v) sin(ξτ ), (13)
For a given v (stationary point, equations ( 10)) and τ we obtain a system of equations with respect to ξ and µ. The value of µ determines the stability boundary and the value of ξ the period of oscillations.

Reaction-diffusion system with delays

Behavior of solutions of system (1), (2) depends on the functions ϕ(v) and ψ(v) which determine production and death of immune cells. We will consider some particular cases below.

Constant ϕ, linear ψ

If we assume ϕ(v) ≡ ϕ 0 (see, e.g. [START_REF] Baker | Computational approaches to parameter estimation and model selection in immunology[END_REF]) and ψ(v) = ψ 0 v reflecting the loss of CTLs on elimination of infection, then for c 0 = 0 we obtain the system

∂v ∂t = D 1 ∂ 2 v ∂x 2 + kv(1 -v) -σvc, ( 14
)
∂c ∂t = D 2 ∂ 2 c ∂x 2 + ϕ 0 c(1 -c) -ψ 0 v µ c. ( 15
)
The initial condition is considered in the form

v(x, t) = 0, c(x, t) = 0.5, 0 < x < L/2; v(x, t) = 0.5, c(x, t) = 0, L/2 < x < L,
where L is the length of the computational interval, -max(τ, µ) ≤ t ≤ 0. We consider the no-flux boundary conditions for both variables.

In the case without delay (µ = 0) system ( 14), ( 15) is similar to the well known system of competition of species in population dynamics. We will consider it also on the whole axis in order to study travelling wave solutions. The corresponding DDE system (without diffusion) has four stationary points, P 0 = (0, 0), P 1 = (1, 0), P 2 = (0, 1), and P 3 = (v * , c * ), where the coordinate of the last point are positive and they can be easily found from the linear algebraic system of equation. The point P 0 is unstable as a solution of the DDE system, points P 1 and P 2 can be stable provided that P 3 is unstable and vice versa. Here and in what follows we say that the point is stable if all its eigenvalues have negative real parts, and that it is unstable if one of the eigenvalues has a positive real part. Consider first the case where the points P 1 and P 2 are stable and the point P 3 is unstable. Then there exist travelling wave solutions of system ( 14), [START_REF] Notkins | Effect of Virus Infections on the Function of the Immune System[END_REF], that is the solution v(x, t) = V (x -st), c(x, t) = C(x -st), where s is the speed of propagation and the functions V , C satisfy the system of equations

D 1 V ′′ + sV ′ + k ′ V (1 -V ) -σV C = 0, ( 16
)
D 2 C ′′ + sC ′ + ϕ 0 c(1 -c) -ψ 0 V (x + sµ)C = 0 ( 17
)
with the conditions

C(-∞) = 1, C(∞) = 0, V (-∞) = 0, V (∞) = 1 (18)
at infinities. Existence of such solution is known for the case without delay (µ = 0) [START_REF] Volpert | Traveling wave solutions of parabolic systems[END_REF] and for a single equation with delay (see [START_REF] Bocharov | Spatiotemporal dynamics of virus infection spreading in tissues[END_REF] and the references therein). We expect that the existence result remains valid for the system of two equations with delay since the methods based on maximum principle remain applicable here.

The speed of the wave can be positive, negative or zero depending on the values of parameters. Figure 4 (left) shows the wave with the zero speed, in other words a stationary solution of system ( 14), [START_REF] Notkins | Effect of Virus Infections on the Function of the Immune System[END_REF]. The speed remains zero in a large interval of variation of time delay µ. Proposition 3.1. Let s µ be the speed of the travelling wave of system ( 14), [START_REF] Notkins | Effect of Virus Infections on the Function of the Immune System[END_REF]. If s 0 > 0 (s 0 < 0), then s µ > 0 (s µ < 0) for all µ ≥ 0.

Proof. For the proof of this assertion it is sufficient to note that s µ is a continuous function of µ. Therefore if s 0 > 0 and s µ0 < 0 for some µ 0 > 0, then s µ = 0 for some 0 < µ < µ 0 . Hence system ( 16), [START_REF] Baker | Computational approaches to parameter estimation and model selection in immunology[END_REF] does not contain delay. The speed of the wave, since it is unique, becomes s µ = s 0 > 0. This contradiction proves the proposition.

Next, consider the case where the point P 3 is stable while the points P 1 and P 2 are unstable. In the case without delay it is known [START_REF] Volpert | Traveling wave solutions of parabolic systems[END_REF] that system [START_REF] Bocharov | Modelling the dynamics of LCMV infection in mice: conventional and exhaustive CTL responses[END_REF], [START_REF] Baker | Computational approaches to parameter estimation and model selection in immunology[END_REF] have solutions with the limits

C(-∞) = 1, C(∞) = c * , V (-∞) = 0, V (∞) = v * ( 19
)
for s ≤ s -, and with the limits

C(-∞) = c * , C(∞) = 0, V (-∞) = v * , V (∞) = 1 ( 20
)
for s ≥ s + . Here s + is some positive and s -some negative values. Numerical simulations show that such solutions also exist for a positive time delay µ (Figure 4, right). The absolute values |s ± | increase as a function of µ.

Linear ϕ and ψ

We now consider linear functions ϕ and ψ: ϕ(v) = ϕ 1 v and ψ(v) = ψ 1 v. The system becomes as follows:

∂v ∂t = D 1 ∂ 2 v ∂x 2 + kv(1 -v) -σvc, ( 21
)
∂c ∂t = D 2 ∂ 2 c ∂x 2 + ϕ 1 v τ c(1 -c) -ψ 1 v µ c. ( 22
)
The corresponding DDE system has the following non-negative stationary points:

v = 0, c ≥ 0; v = 1, c = 0; v * = 1 - σ k ( 1 - ψ 1 ϕ 1 ) , c * = 1 - ψ 1 ϕ 1 .
We consider the same initial conditions as before. For the values of parameters indicated in the caption to Figure 5

D 1 = D 2 = 10 -5 , k = σ = 1, ϕ(v) = 3v, ψ(v) = v; τ = µ = 1 (left), τ = 1, µ = 10, D 1 = D 2 = 10 -6 (middle), τ = 10, µ = 1 (right).
Behavior of solutions depends on the values of time delays. For µ = 10, τ = 1, the stationary point (v * , c * ) remains stable but there are decaying oscillations behind the wave fronts (Figure 5, middle. Decreased diffusion coefficients shorten the space oscillations). If µ = 1, τ = 10, then this point is unstable and we observe spatiotemporal oscillations between the two waves (Figure 5, right). If both time delays are large and equal to each other, µ = 10, τ = 10, then the stationary point becomes again stable.

Bell-shaped ϕ and linear ψ

Finally we will consider a a most general bell-shaped function ϕ(v) describing specifically the induction and down-regulation of immune response and a linear function ψ(v) = ψ 2 v parameterizing the effect of viral population on death of immune cells:

∂v ∂t = D 1 ∂ 2 v ∂x 2 + kv(1 -v) -σvc, ( 23
)
∂c ∂t = D 2 ∂ 2 c ∂x 2 + ϕ(v τ )c(1 -c) -ψ 2 v µ c. ( 24
)
The first three stationary points of the corresponding ODE system are the same as in the case where both functions are linear (Section 3.1): P 0 = (0, 0), P 1 = (1, 0), P 2 = (0, 1). The point P 3 = (v * , c * ) with positive coordinates can be found from the algebraic system of equations:

c = k σ (1 -v), c = 1 - ψ 2 v ϕ(v) .
Note that, its solution can be non-unique. Moreover, stability of stationary points P 1 and P 2 versus P 3 is not mutually exclusive, as it was before. Figures 6 and 7 show different regimes observed in numerical simulations under variation of the parameter σ, specifying the per capita elimination rate of the infection by CTLs and all other parameters fixed. For small values of σ, the virus elimination is not sufficient to stop its propagation, and a single wave with a negative velocity and with the limits

v(-∞) = 0, c(-∞) = 1; v(∞) = 1, c(+∞) = 0 (25)
is observed (Figure 6, left). For a larger value of σ, there are two consecutive waves with negative speed propagating one after another (Figure 6, left). The distance between them grows with time. For σ = 1, the two waves move in the opposite directions (Figure 6, right) Further increase of this parameter leads to the appearance of two consecutive waves with positive speeds (Figure 7, middle). Finally, if it is sufficiently large, there is a single wave with a positive speed (Figure 7, right). It is remarkable that regimes with two consecutive waves propagating with different speeds are specific for the bell-shaped function ϕ(v), they were not observed above (Section 3.2). Time delay influences the wave structure through a possible instability of the steady state with positive coordinates. We take, as example, the result presented in Figure 6 (right) and vary the values of τ and µ. As before, there is a wave with a negative speed and a wave with a positive speed. For τ = 10, µ = 1 they are separated by a growing zone of spatiotemporal oscillations (Figure 8, left). The wave with a negative speed is a solution of the system

D 1 V ′′ + sV ′ + kV (1 -V ) -σV = 0, ( 26
)
D 2 C ′′ + sC ′ + ϕ(V (x + sτ ))C(1 -C) -ψ 2 V (x + sµ)C = 0 ( 27
)
with s < 0 and the limits

V (-∞) = 0, C(-∞) = 1 ; oscillating solution at + ∞. ( 28
)
The wave with the positive speed is a solution of system (26), (27) with s > 0 and oscillating solution at

-∞ ; V (∞) = 1, C(∞) = 0. ( 29
)
For µ = 5 the amplitude of the oscillations decreases, and for µ = 10 they completely disappear.

Discussion

In this work we study spreading of viral infection in the tissues such as lymph node or spleen. Note that the results of our study can be extrapolated to other tissues as the governing equations remain generally valid. The model consists of a reactiondiffusion system of equations with two delays describing the population dynamics of virus and immune cells. The corresponding ODE and DDE systems have different equilibria including a virus free state, which can be stable or unstable, and a disease state where a high virus concentration can be accompanied by a decreased of completely exhausted population of immune cells. The values of parameters and the initial conditions determine the convergence to one of these equilibria.

Virus transport in the tissue plays an important role in the disease development [START_REF] Boulle | HIV Cell-to-Cell Spread Results in Earlier Onset of Viral Gene Expression by Multiple Infections per Cell[END_REF], [START_REF] Mothes | Virus Cell-to-Cell Transmission[END_REF]. Taking it into account by formulation a reaction-diffusion model, we observe the existence of multiple steady-states and possible transition between the different equilibria. Similar to the study based on a single PDE model [START_REF] Bocharov | Spatiotemporal dynamics of virus infection spreading in tissues[END_REF], behavior of the system is determined by the characteristics of antiviral immune response and by the initial viral load.

Proliferation and death of immune cells depends on the infection level. In order to simplify the analysis of the model, we suppose throughout this work a linear dependence of cell death on the concentration of virus. Contrary to this assumption, the dependence of cell proliferation on virus concentration is considered in three different forms. We begin with the simplest case where the proliferation rate does not depend on virus concentration reflecting the model of LCMV infection identified in [START_REF] Baker | Computational approaches to parameter estimation and model selection in immunology[END_REF]. In this case the studied model is similar to the system of competition of species where there are three main types of behavior of solutions:

• (i) The infection can spread in the tissue leading to a complete elimination of immune cells (exhaustion); • (ii) The infection can be completely eliminated (healing); • (iii) There is a stable equilibrium in which virus population co-exists with a population of virus-specific immune cells.

Importantly, in (i) and (ii) the solutions represent a travelling wave propagating either with positive or with negative speed. In the last case (iii), the system converges to the steady-state with an intermediate viral load which can be interpreted as a chronic infection. Next, we consider the case where the proliferation rate of the immune cells is a linear function of the viral load reflecting another mode of immune response regulation as detailed in [START_REF] Baker | Computational approaches to parameter estimation and model selection in immunology[END_REF]. The main difference with the previous case is that the coexisting equilibrium can lose its stability because of the delay in the immune response development. The loss of stability leads to the emergence of spatiotemporal structures. The second time delay in the death term of immune cells is hown that have a stabilizing effect on the steady state.

Finally, the last case is most interesting and biologically realistic. The proliferation rate is represented by a bell shaped function of the virus concentration. It increases for small concentrations and decreases if the concentration is sufficiently high. We observe here two novel regimes not observed before. Infection propagation in the tissue can occur as a superposition of two travelling waves: first wave propagates as a low level infection front followed by a high level infection front with a smaller speed of propagation. Both of the travelling waves can have a positive or a negative speed corresponding to the infection advancement or retreat. As before, these regimes can be accompanied by instabilities and the emergence of spatiotemporal patterns.

Understanding the regularities in the spatio-temporal dynamics of virus infections in host organisms is highly important for improving the clinical practice (e.g. infection fate regulation, biopsy procedures) and in immune inspired engineering (e.g., pattern recognition, classification, control of network systems). There are rather few models with describe the population dynamics of viruses and immune cell using distributed parameter system such as reaction-diffusion equations [START_REF] Bocharov | Spatiotemporal dynamics of virus infection spreading in tissues[END_REF]. In this work, we have analyzed the system of reaction-diffusion PDEs by a combination of analytical and numerical studies. The model is novel in terms of the complexity of the reaction terms describing in a parsimonious way the feedback regulation of activation-, responsiveness tuning-, and apoptosis of immune cells. Our model goes beyond a common view of the virus and immune cell interactions as a predator prey system [START_REF] Nowak | Virus Dynamics[END_REF][START_REF] Bell | Predator-prey equations simulating an immune response[END_REF].

A central problem of infectious disease pathogenesis is related to the mechanisms of T cell exhaustion, associated with various outcomes of infections. Mathematical models assist in revealing a mechanistic view of the infection course and outcomes [START_REF] Bocharov | Understanding Experimental LCMV Infection of Mice: The Role of Mathematical Models[END_REF][START_REF] Nowak | Virus Dynamics[END_REF]. In contrast to previous studies of exhaustion phenomenon, we have shown that in chronic infections the dynamics of virus and immune cells in space and time can be a continuum of different regimes ranging from steady states with a complete or a partial elimination of immune cells on one sisde to non-monotone dynamics in the form of travelling waves and irregular oscillations on the other side. The timedelays reflecting the clonal expansion of cells and the programming for apoptosis appear to be key control parameters. The mapping of the parameter space on the patters of spatio-temporal dynamics being a subject of our study can further assist in resolving via a multi-disciplinary approaches a fundamental problem of how to reverse T cell exhaustion in chronic virus infections [START_REF] De Silva | Resisting immune exhaustion in HIV-1 infection[END_REF][START_REF] Crawford | Molecular and transcriptional basis of CD4 T cell dysfunction during chronic infection[END_REF][START_REF] Speiser | T cell differentiation in chronic infection and cancer: functional adaptation or exhaustion?[END_REF].

Figure 1 .

 1 Figure 1. Schematic representation of the processes underlying the population dynamics of virus infection spread and the antiviral immune response.

  of Parallel, Emergent and Distributed Systems imm-sys3 4 Bocharov, Meyerhans, Bessonov, Trofimchuk, Volpert

Figure 2 .

 2 Figure 2. Schematic representation of the phase plane of system (4), (5). There are 6 stationary points in the left figure. Two of them are stable (see the explanation in the text). The function f (v) becomes negative for large v in the right figure. There are also two stable stationary points, the most left and the most right ones.

  of Parallel, Emergent and Distributed Systems imm-sys3 6 Bocharov, Meyerhans, Bessonov, Trofimchuk, Volpert

  of Parallel, Emergent and Distributed Systems imm-sys3 Parallel, Emergent and Distributed Systems 7

Figure 3 .

 3 Figure 3. Example of periodic oscillations for the values of parameters k = σ = 1, ϕ(v) = 2v, ψ(v) = v, τ = 10, µ = 2 (left). Stability boundary on the (τ, µ) plane for k = σ = 1 (right): Curve 1 is defined by ϕ(v) = 3v, ψ(v) = v, and curve 2 by ϕ(v) = 2v, ψ(v) = v. The stability region is located above the curves.

Figure 4 .

 4 Figure 4. The values of parameters: D 1 = D 2 = 10 -5 , k = σ = 1, ϕ 0 = 1, ψ 0 = 1; τ = µ = 1 (left), σ = 0.5, ψ 0 = 0.5 (right).

  of Parallel, Emergent and Distributed Systems imm-sys3 Parallel, Emergent and Distributed Systems 9

  of Parallel, Emergent and Distributed Systems imm-sys3 10 Bocharov, Meyerhans, Bessonov, Trofimchuk, Volpert

  the stationary point (v * , c * ) with the positive coordinates is stable for sufficiently small time delays τ and µ. In this case, we observe two travelling waves. The first one with the negative speed provides the transition from the point v = 0, c = 0.5 (determined by the initial condition) to the point (v * , c * ), and the second one with the positive speed transition from (v * , c * ) to v = 1, c = 0 (Figure5, left). The interval where the solution converges to the positive equilibrium grows in time.

Figure 5 .

 5 Figure 5. The values of parameters:D 1 = D 2 = 10 -5 , k = σ = 1, ϕ(v) = 3v, ψ(v) = v; τ = µ = 1 (left), τ = 1, µ = 10, D 1 = D 2 = 10 -6 (middle), τ = 10, µ = 1 (right).

  of Parallel, Emergent and Distributed Systems imm-sys3 Parallel, Emergent and Distributed Systems 11

Figure 6 .

 6 Figure 6. The values of parameters: D 1 = D 2 = 10 -5 , k = 1, ψ(v) = 5v; τ = µ = 1; σ = 0.2 (left), σ = 0.5 (middle), σ = 1 (right).

Figure 7 .

 7 Figure 7. The values of parameters: D 1 = D 2 = 10 -5 , k = 1, ψ(v) = 5v; τ = µ = 1; σ = 1.5 (left), σ = 2 (right).

  of Parallel, Emergent and Distributed Systems imm-sys3 12 Bocharov, Meyerhans, Bessonov, Trofimchuk, Volpert

Figure 8 .

 8 Figure 8. The values of parameters: D 1 = D 2 = 10 -5 , k = σ = 1, ψ(v) = 5v; τ = 10, µ = 1 (left); τ = 1, µ = 5 (middle), τ = 10, µ = 10 (right).

  of Parallel, Emergent and Distributed Systems imm-sys3 Parallel, Emergent and Distributed Systems 13

  of Parallel, Emergent and Distributed Systems imm-sys3 14

January 9, 2017 20:38 The International Journal of Parallel, Emergent and Distributed Systems imm-sys3 Parallel, Emergent and Distributed Systems

Acknowledgement

This work was supported by the Russian Science Foundation under Grant [number 15-11-00029].


Andrei Goncharov

Olivier Lopez

Frédéric Anne Amy-Klein

Du Burck

Frédéric Du Burck

Absolute frequency measurements for

Keywords: iodine spectroscopy, hyperfine interaction, optical frequency measurement, femtosecond comb PACS numbers: 06.20.-f, 06.20.F-, 06.30.Bp, 06.30.Ft, 33.20.-t

Introduction

The numerous iodine absorption lines in the 0.5 -0.9 μm spectral range are widely used for laser stabilization in the realization of optical length standards. In particular, seven radiations recommended by the CIPM for the realization of the metre correspond to iodine transitions [START_REF] Quinn | Practical realization of the definition of the metre including recommended radiations of other optical frequency standards[END_REF]. The iodine absorption spectrum is composed of the hyperfine structure of rovibrational transitions between the ground electronic state X (

) and the excited state B ( ). Narrow lines are observed in the spectral range 520-498 nm, near the dissociation limit of the molecule [START_REF] Cheng | Sub-Doppler molecular-iodine transitions near the dissociation limit (523-498 nm)[END_REF]. They are of particular interest as references in compact sources realized with frequency doubled diodepumped solid-state lasers based on Yb

+ Σ g 1 + Π u 0 3 3+
-doped materials which cover this spectral region [START_REF] Wallerand | Absolute frequency measurement of molecular iodine lines at 514.7 nm interrogated by a frequency-doubled Yb-doped fibre laser[END_REF][START_REF] Jacquemet | Blue-green single-frequency laser based on intracavity frequency doubling of a diode-pumped Ytterbium-doped laser[END_REF].

The width of the hyperfine levels in the excited state is limited by their radiative width and by the predissociation effects resulting from the gyroscopic and hyperfine couplings with the dissociative state [START_REF] Broyer | Etude de l'état B de la molécule d'iode par excitation laser : propriétés magnétiques[END_REF][START_REF] Vigué | Fluorescence yield of individual hyperfine components in the B state of I 2 : a test of hyperfine predissociation[END_REF]. The natural predissociation becomes negligible for the high vibrational levels [START_REF] Vigué | Natural hyperfine and magnetic predissociation of the I 2 B state. II -Experiments on natural and hyperfine predissociation[END_REF] and the width of the hyperfine levels are then essentially given by their radiative width, which decreases near the dissociation limit of the molecule [START_REF] Broyer | Etude de l'état B de la molécule d'iode par excitation laser : propriétés magnétiques[END_REF]. Lifetime measurements of rovibronic states near this limit by fluorescence decay have shown that hyperfine transitions with a natural width lower than 10 kHz are expected near 500 nm. In a previous work, we have obtained a half-width at half-maximum (HWHM) experimental linewidth of about 30 kHz for the R(26) 62-0 transition, accessible with an Argon laser at 501.7 nm [START_REF] Du Burck | Recent progress in ultra-high resolution spectroscopy of molecular iodine in Laser Spectroscopy[END_REF] and recently we used such a narrow line detected in a cell for the long term stabilisation of the laser [START_REF] Du Burck | Frequency Measurement of Ar + Laser Stabilized on Narrow Lines of Molecular Iodine at 501.7 nm[END_REF].

+ Π u 0 3 u 1 1 Π
The hyperfine splittings of the R(26) 62-0 transition was first measured by Sorem et al. in the 70's with an uncertainty of about 1 MHz [START_REF] Sorem | Saturation spectroscopy of molecular iodine using the 5017 Å argon laser line[END_REF]. More recently, a detailed study of the iodine hyperfine structure at 514.5 nm and 501.7 nm was reported in [START_REF] Ch | High precision saturation spectroscopy of 127 I 2 with argon lasers at 5145 Å and 5017 Å: I-Main resonances[END_REF]. The frequency differences between the components of the R(26) 62-0 transition were measured with an uncertainty of 5 kHz.

In the present work, we propose a more accurate determination of the hyperfine splittings of the R(26) 62-0 transition. It is based on the direct frequency measurement of an Ar + laser locked successively to every hyperfine components of the transition. They are detected in a continuously pumped cell by saturated absorption spectroscopy at low pressure (0.33 Pa). Thus, narrow resonances with a HWHM of about 50 kHz are used as frequency references. The uncertainty of the frequency measurements is lower than 1 kHz.

Experimental set-up

Some modifications were made to our Ar + /I 2 standard compared with the set-up used for the first absolute frequency measurement of the component a 7 of the R(26) 62-0 transition at 501.7 nm [START_REF] Du Burck | Frequency Measurement of Ar + Laser Stabilized on Narrow Lines of Molecular Iodine at 501.7 nm[END_REF] whereas the system for absolute frequency measurement is essentially the same. We only give here a brief description of the set-up (Fig. 1) and we specify the differences with the previous system.

The single mode Ar + laser is first prestabilized to a Fabry-Perot resonator mode by the Pound-Drever-Hall technique, in order to reduce its frequency jitter. In [START_REF] Du Burck | Frequency Measurement of Ar + Laser Stabilized on Narrow Lines of Molecular Iodine at 501.7 nm[END_REF], it was then locked to a hyperfine component of the R(26) 62-0 transition detected in a low-pressure iodine cell for long term stabilisation. In fact, it was observed that the stability was limited by mechanical vibrations of the Fabry-Perot cavity at a few hertz. To overcome this effect, the laser is now locked to the iodine transition detected in saturated absorption spectroscopy in another cell at a higher pressure (3.3 Pa).

The error signal is obtained by the modulation transfer technique: the pump beam is frequency modulated at 125 kHz and the first harmonic of the saturated absorption signal is detected on the probe beam. In order to improve the signal-to-noise ratio (SNR), the probe beam power is stabilized in a narrow frequency band centred at 125 kHz [START_REF] Du Burck | Narrow-band correction of the residual amplitude modulation in FM spectroscopy[END_REF][START_REF] Du Burck | Frequency Measurement of Ar + Laser Stabilized on Narrow Lines of Molecular Iodine at 501.7 nm[END_REF].

The Ar + laser, prestabilized to the Fabry-Perot resonator mode and the high-pressure iodine cell is then locked to the same hyperfine component of the R(26) 62-0 transition detected in the low-pressure iodine cell. The resonance is detected in saturated absorption spectroscopy with the technique of frequency modulation (FM) spectroscopy [START_REF] Bjorklund | Frequency-modulation spectroscopy: a new method for measuring weak absorptions and dispersions[END_REF][START_REF] Hall | Optical heterodyne saturation spectroscopy[END_REF]. The probe beam is frequencymodulated at with an acousto-optic modulator (AOM). The residual amplitude modulation generated by the AOM is cancelled by the narrow-band control of the beam intensity at the modulation frequency [START_REF] Du Burck | Narrow-band correction of the residual amplitude modulation in FM spectroscopy[END_REF]. The phase sensitive detection of the probe beam at 2.5 MHz gives an error signal with a dispersion shape. An additional amplitude modulation of the saturation beam at and a second phase sensitive detection are used to eliminate the residual offset of the error signal due to the Doppler background. The length of the low pressure cell is 4 m, its diameter is 15 cm. The beam (diameter 6 mm) passes two times through it to reach an 8 m interaction length with the iodine vapour. This cell is pumped during the experiment to minimize buffer gases and impurities effects. Indeed, the width of the studied transitions at 501.7 nm is known to be very sensitive to the impurities present in the cell [START_REF] Flory | Perturbation de l'iode moléculaire en phase vapeur par la présence de gaz étrangers, conséquence pour la stabilisation du laser He[END_REF]. The iodine pressure is controlled by thermostabilization of a cold finger. The narrow lines detected with this cell in previous works using Raman and Rayleigh spectroscopy [START_REF] Du Burck | Resolving power and sensitivity in modulation transfer stimulated resonant Raman spectroscopy[END_REF][START_REF] Du Burck | High-frequency modulation transfer technique for ultra-high resolution spectroscopy of I 2[END_REF] show that the upper limit of impurity pressure in the cell is less than 0.04 Pa (0.3 mTorr).

The frequency of the Ar + laser locked to the studied hyperfine component is compared to a mode of a wide bandwidth optical comb generated by a Kerr-lens mode-locked Ti:Sa laser spectrally enlarged by a photonic crystal fibre. One detects the beat note signal between a small amount of the laser power (about 10 mW) and the optical comb. It is then electronically shifted by the carrier envelop offset frequency, detected with the usual f-2f set-up, and is measured with a frequency counter (for details, see [START_REF] Du Burck | Frequency Measurement of Ar + Laser Stabilized on Narrow Lines of Molecular Iodine at 501.7 nm[END_REF]). A local oscillator locked to a 1 GHz reference drives the frequency counter and the repetition rate of the femtosecond laser. The reference signal with a high level of stability and accuracy (cryogenic oscillator/H-maser/Cs-fountain) [START_REF] Chambon | Design and realization of a flywheel oscillator for advanced time and frequency metrology[END_REF] is received from LNE-SYRTE at Observatoire de Paris through a 43 km long optical fibre, as an amplitude modulation of a 1.5 μm optical carrier. The whole frequency measurement set-up has a resolution of a few 10 -14 for 1 s integration time [START_REF] Daussy | Long-distance frequency dissemination with a resolution of 10 -17[END_REF].

Measurement of the hyperfine structure of the R(26)62-0 transition

The frequency standard at 501.7 nm shows a stability better than 10 -12 for 1 s integration time, limited by the signal to noise ratio of the detected line in the low pressure cell. Several factors can affect the frequency reproducibility. The frequency shifts with the iodine pressure, the beam power and the modulation index were investigated in [START_REF] Du Burck | Frequency Measurement of Ar + Laser Stabilized on Narrow Lines of Molecular Iodine at 501.7 nm[END_REF].

-The iodine pressure in the cell is controlled by the temperature of a cold finger. The temperature is -22.5°C corresponding to an iodine pressure of 0.33 Pa. The HWHM linewidth of the resonance is then 45 kHz (Fig. 2). The pressure shift coefficient is less than . At -22.5°C, the sensitivity of the iodine pressure with temperature is 0.04 Pa/K, which gives a very low sensitivity of the measured frequency with the temperature of about -150 Hz/K. kHz/Pa 8 . 3 --The probe beam power is 0.5 mW. Its decrease by a factor 2 leads to a shift of about 1 kHz.

That corresponds to a sensitivity of the frequency with the power of 4 kHz/mW.

-The modulation index is close to unity which is the optimum value for the signal detection.

The modification of the modulation index does not lead to any appreciable shift.

Considering an uncertainty of the cold finger temperature of 1 K, the frequency uncertainty associated to the iodine pressure is 150 Hz. The uncertainty of the beam power is estimated to be In a first step, several measurements of the absolute frequency of component a 7 of the R(26) 62-0 transition were repeated during two successive days. It is seen in Fig. 3 that the deviations between measurements operated from one day to the other does not exceed the deviations observed during measurements operated on the same day ( Hz 250 ≤

). The mean value of the 8 absolute frequency measurements of component a 7 is found to be 597 366 498 654.47 kHz.

The standard deviation is 160 Hz, a smaller value than the combined standard uncertainty of 250 Hz estimated above.

In a second step, the laser was locked successively on the fifteen hyperfine components of the R(26) 62-0 transition and their absolute frequencies were measured. The absolute frequency of component a 1 is found to be 597 366 101 477.9 kHz. The positions of the other lines, referred to component a 1 are reported in Table 1.

The hyperfine splittings of the R(26) 62-0 transition deduced from our absolute frequency measurements can be compared to data reported in [START_REF] Ch | High precision saturation spectroscopy of 127 I 2 with argon lasers at 5145 Å and 5017 Å: I-Main resonances[END_REF]. In this work, the frequency differences between hyperfine components were obtained from the differences of the beat frequency of two Ar + lasers, one of them locked to a reference component, and the other one locked successively to the hyperfine components of the transition. We calculated the difference between the splittings deduced in the present work from absolute frequency measurements and those given in [START_REF] Ch | High precision saturation spectroscopy of 127 I 2 with argon lasers at 5145 Å and 5017 Å: I-Main resonances[END_REF]. The global frequency shift between both sets of data is compensated by taking the sum of the deviations equal to 0. The resulting values are given in the last column of Table 1. The maximum deviation between both sets of data is 12.2 kHz and the standard deviation is 6.8 kHz. This latter value is of the same order as the uncertainty of the frequency differences (5 kHz -1 σ) estimated in [START_REF] Ch | High precision saturation spectroscopy of 127 I 2 with argon lasers at 5145 Å and 5017 Å: I-Main resonances[END_REF] from the statistical reproducibility.

Conclusion

We have performed for the first time, the absolute frequency measurement of the complete hyperfine structure of the R(26) 62-0 transition in molecular iodine at 501.7 nm with an optical clockwork based on a femtosecond laser frequency comb generator. Our Ar + laser was locked to narrow resonances of about 50 kHz linewidth, detected in a low pressure continuously pumped cell with a low impurity pressure. The uncertainty of the frequency measurements is Hz 250 = σ . These measurements give a better accuracy on the hyperfine splittings of the R(26) 62-0 transition. These data can be combined with the measurement of hyperfine splittings of the lower

Table captions

Tab. 1: Absolute frequencies of the hyperfine components of the R(26) 62-0 transition in iodine and deviations between the frequencies deduced from [START_REF] Ch | High precision saturation spectroscopy of 127 I 2 with argon lasers at 5145 Å and 5017 Å: I-Main resonances[END_REF] and our frequencies, taking the sum of the deviations equal to 0.

Figure captions

  10 %, corresponding to a frequency uncertainty of 200 Hz. The frequency measurement chain does not introduce any frequency error. The combined standard uncertainty (1σ) is then

Fig. 1 :

 1 Fig. 1: Ar + /I 2 frequency standard at 501.7 nm.

Fig. 2 :

 2 Fig. 2:Error signal for the long term stabilization of the Ar + laser on the hyperfine component a 7

Fig. 3 :

 3 Fig. 3: Frequency measurements of the hyperfine component a 7 of the R(26) 62-0 transition. The

Fig. 3

 3 Fig.1

Table 1 component

 1

		f(a n ) -597 366 101 477.9 kHz	Frequencies deduced
		(kHz)	from [11] -ours
			(kHz)
	a 1	0	6.5
	a 2	166 253.5	5.0
	a 3	265 934.0	1.6
	a 4	295 233.9	-4.4
	a 5	317 938.2	-1.6
	a 6	364 499.5	-11.0
	a 7	397 176.6	11.0
	a 8	474 039.1	-0.6
	a 9	523 548.2	-6.7
	a 10	560 031.4	-2.9
	a 11	637 054.1	-5.6
	a 12	660 262.7	3,8
	a 13	739 965.5	1.0
	a 14	769 747.2	-8.7
	a 15	842 244.3	12.2

Acknowledgments

The authors wish to thank Ch. J. Bordé for his interest in this work. This work was supported by GDRE project "Lasers et techniques de l'information" and by RFBR grant #05-02-19645.


Marc Shapiro

Conflict-Free Replicated Data Types (basic entry)

Keywords: Commutative Replicated Data Types (CmRDTs), Convergent Replicated Data Types (CvRDTs), Replicated

come

Historical background

There is precursor work on specific CRDTs, before the concept was formally identified as an independent abstraction. Johnson and Thomas [START_REF] Johnson | The maintenance of duplicate databases[END_REF] proposed the so-called Last-Writer-Wins (LWW) or Greatest-Timestamp-Wins approach for a replicated register, 1 i.e., an untyped memory that an update completely overwrites. Wuu and Bernstein [START_REF] Wuu | Efficient solutions to the replicated log and dictionary problems[END_REF] studied more complex data types, the log and dictionary (a.k.a. map or Key-Value Store). The whole area of Operational Transformation (OT) studied replicated strings or sequences, intended for concurrent editing applications [START_REF] Sun | Operational transformation in real-time group editors: issues, algorithms, and achievements[END_REF]. Baquero and Moura [START_REF] Baquero | Using structural characteristics for autonomous operation[END_REF] identified some convergence conditions for data types used in mobile computing. The Dynamo system is based on a multi-value register construct [START_REF] Decandia | Dynamo: Amazon's highly available key-value store[END_REF]. Related topics include replicated file systems and version control systems.

The concept of CRDTs was identified by Preguiça et al. [START_REF] Preguiça | A commutative replicated data type for cooperative editing[END_REF], formalised by Shapiro et al. [START_REF] Shapiro | Conflict-free replicated data types[END_REF] and Shapiro et al. [START_REF] Shapiro | Convergent and commutative replicated data types[END_REF], and studied systematically in Shapiro et al. [START_REF] Shapiro | A comprehensive study of Convergent and Commutative Replicated Data Types[END_REF]. A similar concept, called Replicated Abstract Data Types (RADTs), was proposed independently by Roh et al. [START_REF] Roh | Replicated Abstract Data Types: Building blocks for collaborative applications[END_REF]. These works consider symmetric replicas, in which concurrent updates must be commutative and associative. In related work, Burckhardt et al. [START_REF] Burckhardt | Cloud types for eventual consistency[END_REF] consider so-called Cloud Types with asymmetric main and secondary branches, thus relaxing the commutativity requirement.

The distinction between operation-and state-based CRDTs was established by Shapiro et al. [START_REF] Shapiro | Conflict-free replicated data types[END_REF]. Burckhardt et al. [START_REF] Burckhardt | Replicated data types: Specification, verification, optimality[END_REF] established lower-bound and optimality results for some representative state-based CRDTs. Delta-CRDTs were proposed to decrease the footprint of state-based CRDTs while keeping most of their advantages [START_REF] Almeida | Efficient state-based CRDTs by deltamutation[END_REF]. Pure operation-based CRDTs leverage causal-order delivery to streamline the design and implementation of operation-based CRDTs [START_REF] Baquero | Making operation-based CRDTs operationbased[END_REF].

Foundations

Encapsulating replication and concurrency

In a distributed system, shared data is often replicated to improve the availability and latency of reads. However, requiring strong consistency between replicas will actually degrade the availability and latency of writes. According to the CAP Theorem, to improve write availability and latency requires to relax the consistency requirement: a replica should accept updates without synchronising with other replicas, and propagate them in the background.

If multiple replicas accept updates (a so-called "multi-master" system), inevitably, there will be concurrent updates to separate replicas. Managing and reconciling conflicting concurrent updates, in order to ensure Eventual Consistency, is a major issue of such systems.

CRDTs were invented to resolve this issue, by encapsulating a familiar object abstraction with a mathematically-sound conflict resolution protocol.

CRDT behaviour

A CRDT supports the interface of the corresponding abstraction. Thus, a register CRDT will support mutation methods such as read and write methods; a counter supports increment, decrement and value methods; a set methods to add, remove and query elements, and so on.

A number of CRDT types have been proposed in the literature. The most basic ones are the LWW Register [START_REF] Johnson | The maintenance of duplicate databases[END_REF] and the Multi-Value Register [START_REF] Decandia | Dynamo: Amazon's highly available key-value store[END_REF]. A widely-studied CRDT is the sequence or list, used in particular for cooperative editing [START_REF] Preguiça | A commutative replicated data type for cooperative editing[END_REF][START_REF] Roh | Replicated Abstract Data Types: Building blocks for collaborative applications[END_REF]. Other common CRDTs include counters [START_REF] Shapiro | A comprehensive study of Convergent and Commutative Replicated Data Types[END_REF], sets [START_REF] Shapiro | A comprehensive study of Convergent and Commutative Replicated Data Types[END_REF] and maps [START_REF] Roh | Replicated Abstract Data Types: Building blocks for collaborative applications[END_REF].

Consider for instance a set data type, supporting operations to add and remove elements, ignoring duplicates. An archetypical CRDT set is the so-called Observed-Remove Set (OR-Set). In any sequential execution, it behaves exactly like a sequential set. Concurrently adding and removing different elements e and f , or adding the same element e twice, or removing the same element e twice, commute per the sequential specification. However, to ensure commutativity of two updates that concurrently add and remove the same element e, the OR-Set makes the "add win," i.e., any replica that observes both operations concludes that e is a member of the set. To do this, the implementation of the remove operation effectively cancels out those add operations that it previously observed, and only those. We return to this example later in this entry.

Implementation approaches and requirements

A CRDT is typically designed to behave like its sequential counterpart in any sequential execution. Furthermore, a CRDT is replicated, supports concurrent updates for availability, and encapsulates some strategy to merge concurrent updates and ensure that its replicas eventually converge. One such strategy is the "Last-Writer-Wins" approach [START_REF] Johnson | The maintenance of duplicate databases[END_REF] that uses timestamps to totally order updates and discard all but the highest-timestamped one. Another is to record concurrent updates side-by-side, so that the application can deal with them later, as in the Multi-Value Register of Dynamo [START_REF] Decandia | Dynamo: Amazon's highly available key-value store[END_REF] and in many filesystems or version control systems.

The literature distinguishes two implementation strategies for CRDTs. In the state-based approach, a mutation method changes only the state of the origin replica. Periodically, a replica sends its full state to some other. The receiver merges the received state into its own. A state-based CRDT manages its state space as a join-semilattice, where every mutator method is an inflation, and the merge method computes the join (a.k.a. least-upper-bound) of the states to be merged [START_REF] Baquero | Using structural characteristics for autonomous operation[END_REF][START_REF] Shapiro | Conflict-free replicated data types[END_REF]. Semi-lattice join is associative, commutative and idempotent. The first two properties ensure that all replicas converge deterministically to the same outcome. The latter ensures that the system tolerates duplicated merges. As long as replicas communicate their state sufficiently often, and the communication graph is connected, replicas eventually converge and each object's history is causally consistent.

The operation-based approach consists of sending updates rather than states. A mutator method consists of two steps. The generator step reads the state of the origin replica and generates an effector, a state transformation that is sent and eventually applied to all replicas in the second step [START_REF] Gotsman | Cause I'm strong enough: Reasoning about consistency choices in distributed systems[END_REF][START_REF] Li | Making geo-replicated systems fast as possible, consistent when necessary[END_REF][START_REF] Shapiro | Conflict-free replicated data types[END_REF]]. 1 Concurrent effectors must commute with one another since they may be received in any order. Associativity is not required but, if available, enables batching multiple effectors into a single one. The operation-based approach requires that the underlying communication layer deliver updates to the object in causal order, and never deliver the same (non-idempotent) update twice.

The state-based approach is generally considered less efficient (state may be very large) but more elegant and simpler to understand. It makes very few assumptions about the underlying network; for instance, the number and identity of replicas may be unknown and variable. Conversely, the operation-based approach appears more efficient but is more complex to implement and requires a more elaborate communication layer.

While CRDTs replicas are guaranteed to eventually converge, this may be insufficient for application correctness. Many applications also require Causal Consistency to avoid ordering anomalies across objects. Furthermore, maintaining the integrity of structural invariants may require synchronisation to disallow certain concurrent updates [START_REF] Balegas | Putting consistency back into eventual consistency[END_REF][START_REF] Gotsman | Cause I'm strong enough: Reasoning about consistency choices in distributed systems[END_REF] (see also Multi Datacenter Consistency).

Example: OR-Set

The following pseudocode, in the style of Shapiro et al. [START_REF] Shapiro | Conflict-free replicated data types[END_REF], illustrates a state-based implementation of an OR-Set [START_REF] Shapiro | A comprehensive study of Convergent and Commutative Replicated Data Types[END_REF]. The local variables of a replica are a set E of added elements, and a set T of removed elements or tombstones. Adding an element e puts it into E along with a unique tag. The tag remains internal to the implementation and is not visible through the interface. Removing an element e moves all pairs of the form (e, ) from S into T . An element e is contained in the set if there exists a pair of the form (e, ) in E. Merging two states retains the element pairs that are contained in both states, and makes tombstones of element pairs that are tombstones in either state. Unfortunately the memory usage of this specification grows, without bound, with every add and remove operation. However, observe that adding an element pair necessarily happens-before removing the same pair. Leveraging this observation, Bieniusa et al. [6] propose an implementation whose size is bounded by the number of currently-contained elements; since the state-based approach does not assume any particular delivery order, it is somewhat complex. As the operation-based approach already assumes causal-order delivery, avoiding tombstones is straightforward, as shown next [START_REF] Shapiro | Convergent and commutative replicated data types[END_REF]. A replica maintains a set of contained element-pairs E. Adding an element e creates the corresponding pair, and removing an element e simply removes all pairs of the form (e, ) observed at the origin replica. Several implementations of CRDTs have been reported, in languages such as C++, Clojure, Erlang, Go, Java, Python, Ruby, and Scala.

The Riak NoSQL database, as of Version 2.0, implements a number of replicated data types, including flags, registers, counters, sets and maps [START_REF] Basho | Data types, version 2.1.1[END_REF]. Bet365, a large online betting company, which manages 2.5 million simultaneous users with Riak OR-Sets. League of Legends, an online multiplayer game, implements online chat for 70 million users with Riak sets. TomTom extend Riak's CRDTs to share navigation data between a user's different devices. SoundCloud uses a Go implementation on top of the Redis database to store time-series information.

-

  -State-based OR-Set specification, with tombstones variables set E, set T --E: elements; T : tombstones --sets of pairs { (element e, unique-tag n), . . . } initial ∅, ∅ query contains (element e) : boolean b let b = (∃n : (e, n) ∈ E) update add (element e) let n = unique() -unique() returns a unique tag E := E ∪ {(e, n)} -e + unique tag update remove (element e) let R = {(e, n)|∃n : (e, n) ∈ E} --Collect all unique pairs containing e E := E \ R T := T ∪ R --Make pairs observed at origin into tombstone merge (B) E := (E \ B.T ) ∪ (B.E \ T ) T := T ∪ B.T

-

  -Operation-based Observed-Remove Set, without tombstones variables set E --set of pairs { (element e, unique-tag u), . . . } initial ∅ query contains (element e) : boolean b let b = (∃u : (e, u) ∈ E) update add (element e) generator (e) let u = unique() -unique() returns a unique value effector (e, u) E := E ∪ {(e, u)} -e + unique tag update remove (element e) generator (e) let R = {(e, u)|∃u : (e, u) ∈ E} --Generator: Collect all unique pairs containing e effector (R) E := E \ R --Effector: remove pairs observed at source 5 Usage

This vocabulary is not standardised. Other names for the generator phase are upstream, prepare, or prepare-update. Alternative names for the effector phase are downstream, effect, effect-update or shadow operation.