A mathematical model for tilapia lake virus transmission with waning immunity

The goal of this paper is to investigate the influence of the waning immunity on the dynamics of Tilapia Lake Virus (TiLV) transmission in wild and farmed tilapia within freshwater. We formulate a model for which susceptible individuals can contract the disease in two ways: (i) direct mode caused by contact with infected individuals; (ii) indirect mode due to the presence of pathogenic agents in the water. We obtain an age-structured model which combines both age since infection and age since recovery. We derive an explicit formula for the reproductive number and show that the disease-free equilibrium is locally asymptotically stable when, . We discuss on the form of the waning immunity parameter and show numerically that a Hopf bifurcation may occur for suitable immunity parameter values, which means that there is a periodic solution around the endemic equilibrium when, .


Introduction
The immune system can be categorized into specific (adaptive) immunity with immunological memory and non-specific innate immunity [14,23,34]. The specific immune system involves the recognition of specific antigens on a pathogen, thereby providing protection against that pathogen primarily in the form of specific antibodies. In contrast, the non-specific immune system provides an array of protective mechanisms that are inherently available and provide immediate and permanent protection against a wide variety of pathogen. Therefore, the adaptive immune response is activated by the innate immune system in order to allow a specific response [31].
Further, the innate parameters are at the forefront of host immune defense and are a crucial factor in disease resistance. The adaptive response is commonly delayed but is essential for long-lasting immunity and a key factor in successful vaccination. Both innate and adaptive immunities complement each other in the host's attempt to prevent infection [14]. A key feature of the adaptive immune system is immunological memory. Repeat infections by the same virus or bacteria induce a strong and specific response that usually stops the infection and has less reliance on the innate system. Vaccination against infection is possible due to this immune memory. The first adaptive response against an infection, called the primary response, often takes days to mature. In contrast, a memory response develops within hours of infection. Memory is maintained by a subset of B and T lymphocytes called memory cells, which can potentially survive for years in the body. Memory cells remain ready to respond rapidly and efficiently to a subsequent encounter with a pathogen, giving rise to stronger and faster so-called secondary responses [20]. Most animals are susceptible to the initial infection, but they differ in their ability to limit the infection or destroy the pathogen [29].
Tilapia is considered as the second largest farmed fish species after carp globally. With such expansion, there is the risk of emergence of infectious diseases and tilapia lake virus infection has been shown to cause high mortalities ( 20-90%) over the last 4-5 years [4,12,25]. The identified routes of transmission are: horizontally through cohabitation [13], and vertically, that is from infected fish to their offspring [11,35]. Recently Tattiyapong et al. in [33], showed that the tilapia can develop protective immunity including a humoral response following exposure to TiLV. In this paper, we devote our attention to the modelling of an age-structured model for TiLV transmission that incorporates both horizontal direct (contact), indirect (through water) transmission routes, and waning immunity. To the best of our knowledge, only two mathematical papers have been devoted to understand the dynamics of the TiLV transmission. Their authors use various mathematical approaches: ordinary differential equation [36] and partial differential equation [19]. The model developed here can be viewed as an extension of the one in [19] in the following sense: • According to [9,18], the death rate in tilapia fisheries was reported to be densitydependent. So, to be more realistic, we consider two types of mortality rates: the density non-dependent and the density-dependent mortality. • Though the dynamics of immune status among host individuals plays a crucial role in the spread of infectious disease, Tattiyapong et al. [33] showed that TiLV infection provides recovered individuals with a short or relatively long immunity against re-infection and they lose their immunity through time, allowing them to become susceptible again. This means that it is natural for us to include the effects of immunity into our model in order to better represent the actual dynamics of TiLV epidemic spread and predict future outbreaks.
Furthermore, several papers have investigated diseases transmission with age of infection [3,7,8,19,21,24,37]. But, only few have addressed models with both time since infection and waning immunity. S. Bhattacharya and F.R. Adler [5] considered a SIRS model with temporary immunity when the rate of loss of immunity can depend on the time since recovery from disease. They determined the conditions under which the endemic steady state becomes unstable and found that rapid immunity loss is necessary to produce oscillations. In [17], H. Gulbudak et al. developed an immuno-epidemiological model that links the within-host dynamics to between-host circulation of a vector-borne disease. They derived evolutionary optimization principles for both pathogen and host and obtained by the invasion analysis that the R 0 maximization principle holds for the vector-borne pathogen. For the host, they find that increasing the vector inoculum size increases the pathogen R 0 , but can either increase or decrease the pathogen virulence, suggesting that vector inoculum size can contribute to the virulence of vector-borne diseases in distinct ways.
More recently, R. Djidjou-Demasse et al. [30] considered a human-vector malaria transmission model structured by age, time since infection and waning immunity. They proved the existence of equilibria and obtained a necessary and sufficient condition that implies the bifurcation of the endemic equilibria. They also proved that by neglecting the age dependence of the human population, and under some certain conditions, they may be a backward or a forward bifurcation.
The present article considers an age-structured model with loss of immunity, which is considered to be one of the sources causing recurrence of infectious disease dynamics observed in many epidemics. Every recovered individual losses progressively its immunity, implying its direct reversion from the recovered class to the susceptible class at a rate that depends upon time since recovery. The obtained system is a SIVRS model, including two hyperbolic first-order partial differential equations coupled with two ordinary differential equations. Using the integrated semi-group theory [22,24], we prove some well-posedness results. We obtain an explicit formula for the basic reproduction number R 0 , we show the existence of a unique disease-free equilibrium and obtain its local stability. Finally, we explore the role of waning immunity on the epidemic evolution.
The organization of this paper is as follows. In Section 2, we formulate the model and prove its mathematical well-posedness. Section 3 is devoted to the computation of the basic reproduction number, the disease-free equilibrium and its stability. In Section 4, we perform some numerical simulations. In Section 5, a brief discussion is given to end this work.

Model formulation
Let S(t) be the density of susceptible fish at time t, i(a, t) the density of fish infected at time t with respect to age of infection a and r(τ , t) be the density of recovered fish at time t with respect to the time since recovery τ . We denote by V(t), the pathogen concentration in a water source at time t. The birth rate in the population is denoted by b. The constant σ ind denotes the transmission rate parameter for water-to-fish contact (horizontal indirect transmission). The density non-dependent and dependent mortality rates are denoted by μ 0 and μ 1 respectively. During an infection, fish can either die from infection at a rate μ i (a), or recover at a time since infection a with a rate δ(a). A recovered fish loses its immunity at a rate γ (τ ). The proportion of offspring born in the infected class is denoted by π(a) ∈ [0, 1], and we assume that offspring from the recovered class are born susceptible.
Then the terms b β(a)i(a, t) da, we assume that a fraction p ∈ [0, 1] of this flux, is directly ingested by fish at time t by direct contact (from fish to fish) with a rate σ dir and the remaining 1−p is shed in water. The flux of newly infected fish corresponds to the boundary condition for i at age a = 0, while the flux of newly recovered fish corresponds to the boundary condition for r at τ = 0. The mortality rate of the virus is denoted by μ V . The total number of infected individuals is Then the total size of population at time t is given by N(t) = S(t) + I(t) + R(t). With all the above notations, the model reads as follows where with boundary conditions and the following initial conditions Assumption 2.1: We make the following assumptions and γ (τ ) > γ 0 , for a.e τ ∈ (0, ∞) and for some γ 0 > 0; (iii) The functions β(.), π i (.) are positive, bounded and uniformly continuous on (0, ∞).

Mathematical well-posedness
In this section, we prove that the system (1)-(4) is well posed by using an integrated semigroup approach [21,24]. Before going further, we introduce the space X 1 defined by Let A i : D(A i ) ⊂ X 1 → X 1 and A r : D(A r ) ⊂ X 1 → X 1 be two linear operators on X 1 defined by Next, let X be the space defined by It is clear that X endowed with the norm · X is a Banach space. We denote by X + the positive cone of X, that is Then We consider F : D(A) → X the nonlinear map defined by Set X + 0 = X 0 ∩ X + . Then, the system (1)-(4) rewrites as the following abstract Cauchy problem where In general, the differential Equation (5) may not have a strong solution. Thus, we solve it in integrated form: A continuous solution to (5) is called an integral solution to (6). The next result states the well posedness of (1)-(4). Proof: We proceed as in [24]. First, it is easy to see that the function F is Lipschitz continuous on bounded sets. Next, we prove that the operator (A, D(A)) is a Hille-Yosida operator and Solving this latter system of equation lead us to Integrating the equations forφ and forψ with respect to a and τ , respectively, and adding all obtained equations yields Hence A is Hille-Yosida. Moreover, if we assume that f ∈ X + , then by (7), we obtaiñ ϕ ≥ 0,φ ≥ 0,ω ≥ 0 andω ≥ 0. That is g = (φ, (0,φ), (0,ψ),ω) ∈ X + . Thus (λ − A) −1 maps X + into X + . Next we show that for all t ≥ 0 and v ∈ X + 0 , where dist(v, with α being chosen such that α > C. Then F(v) ∈ X + and for any positive and sufficiently small h, we have α − αhv ∈ X + and h F(v) ∈ X + . So lim In what follows, we prove that the solution of (1)-(4) is bounded.
Moreover, the upper bounds are uniform, Proof: It is clear that N(t) satisfies the differential inequality Hence it follows that from which we deduce the first inequality of (9). In addition, we obtain Moreover from the third equation of (1), we deduce that Thus from which we deduce the second inequality of (9). Hence This completes the proof.

Remark 2.1:
The Theorem 2.2 shows that the model has a unique continuous solution in a positive cone X + 0 , moreover, the set D defined as follows is positively invariant for system (1).
The next section concerns the basic reproduction number.

Basic reproduction number and disease-free equilibrium
We will consider the following functions and where η 1 (a) = δ(a) + μ i (a) + b and η 2 (τ ) = γ (τ ) + b. Then, it is clear that where R Hdir where F(a) is given by (12).
Moreover the third equation of (20) gives with F(a) given by (12).
Taking (21) and (22) into account, the third equation of (20) leads us tõ That is Define a function Then, H(0) = R 0 . Moreover, it is easy to see that, H is a continuously differentiable satisfying Therefore H is a decreasing function. Hence, any real solution of Equation (23) is negative if R 0 < 1, and positive if R 0 > 1. Thus, if R 0 > 1 the infection-free equilibrium is unstable. Next, we show that Equation (23) has no complex solutions with nonnegative real part if R 0 < 1. To do this, let λ = x + iy, with x, y ∈ R be a solution of Equation (23). Now define and So We argue by contradiction by assuming that x ≥ 0. Then, Mean duration of the laying being affected day 13 [13,19,33] That is |H(λ)| < 1, which is a contradiction. Thus, x < 0 and every solution of (23) has a negative real part. Therefore, if R 0 < 1, the disease-free equilibrium E 0 is locally asymptotically stable and is unstable if R 0 > 1. This completes the proof.

Numerical experiments
In this section, we present some numerical simulations. The list of the parameters of our model as well as their values are summarized in Table 1.
Values of γ (.): to choose a suitable form of γ (τ ), we assume that the duration of the immunity loss decays linearly with time since recovery. Furthermore, according to [33], tilapia developed an antibody response as early as 7 days post TiLV challenge (dpc), peaked at 15 dpc, showed a gradual decline up until about 42 dpc, but persisted in some fish up until day 110 dpc. Hence we estimate that the mean time to loss immunity isT ≈ 90 days. Let (τ ) = e − τ 0 γ (s) ds denote the probability that an immune individual remains immune at time τ after recovery. Then,T Next, we choose γ (τ ) as an increasing function of τ in the form: In what follows, ς will be interpreted as a bifurcation parameter. Since the mean time to loss immunity is estimated toT ≈ 90 days, we get ς = 180.
In the same way, we adopt the following form of δ(a): where ε is set to ε = 0.338 such that the mean duration of the infectious period, that is ∞ 0 e − a 0 δ(s) ds da is around 13 days [33]. Value of parameter μ 1 : The value of the density-dependent mortality μ 1 can vary depending on the carrying capacity. The value of μ 1 was fixed to 8.22 × 10 −8 day −1 . In that case, the total size of fish converges towards the carrying capacity b−μ 0 μ 1 = 10 5 fish at the disease-free equilibrium.
Values of parameters σ dir and σ ind : the Horizontal direct transmission rate σ dir was estimated to 6.8 × 10 −9 copies −1 and the Horizontal indirect transmission rate σ ind to 3.6 × 10 −11 copies −1 .day −1 . These values should be estimated by fitting the model with the experimental data. But unfortunately, we do not have data to estimate these parameters. Here, the values of σ dir and σ ind were calibrated in order to have a reproductive number R Hdir 0 or R Hind 0 very close to that estimated in [19,36], that was R 0 = 2.60 ± 0.16. The following forms of parameters were retrieved from [19]. Values of β(.): Values of π(.): we consider π(a) on the form π(a) = π 0 (π 1 − a)(a − (π 1 + π 2 )) if a ∈ [π 1 , π 1 + π 2 ], 0 o t h e r w i s e ,

Effects of waning on the disease transmission
We present the numerical simulation results and show numerically the existence of endemic equilibria and the persistence of the disease when, R 0 > 1. The simulations are performed with a set of parameters which are described in Table 1. Let We investigate the effect of waning on the spread of the disease. We consider indeed two scenarios: The first one when the waning parameter γ is given in (28) with the bifurcation  Table 1. The waning parameter γ is given in (28) with the bifurcation parameter ς that varies. We show to simplify the graphical representations the quantities

Discussion and conclusion
This work concerns a model of TiLV occurring in tilapia's populations. Because understanding the role of the adaptive immune response following exposure of tilapia to TiLV is a critical step in the development of a vaccine [33], our model incorporates the rate at which recovered individuals lose their immunity. We have established a mathematical wellposedness result obtained by mean of integrated semi-group theory, have computed the basic reproduction number and proved the stability of the disease-free equilibrium. Finally, we have performed some numerical simulations to prove the existence of the endemic equilibrium, the persistence of the disease and to illustrate the role of waning on the epidemic evolution.
The form of R 0 given in (15) has a biological meaning; indeed R 0 depicts the expected number of secondary infections resulting from a single primary infection into an otherwise susceptible population. The term R Hdir 0 is the average number of secondary infections produced by one infective individual during its infectious period by horizontal direct transmission while R Hind 0 represents the average number of secondary infections produced by one infective individual during its infectious period by horizontal indirect transmission. R V 0 is the average number of secondary infections produced by one infective individual during its infectious period by Vertical transmission. Our study showed the possibility of several outcomes depending on the basic reproduction number R 0 , which led us to find mathematically the existence and the local stability of the disease-free equilibrium when, R 0 < 1. Then, using numerical simulations, we found the existence of a unique endemic equilibrium for the system (1)-(4), and the persistence of the disease in the population, when R 0 > 1.
In order to have a better insight of the influence of waning on the dynamics of the model, we performed numerical simulation with waning parameter γ (.) that varies with respect to a bifurcation parameter ς. The results of these simulations in both cases suggest that endemic equilibrium may lose its stability via a Hopf bifurcation, thus giving rise to stable periodic solutions. Biologically, this means that when the temporary immunity period is within a certain range, there will be periodic outbreaks of epidemic, and the disease will not be eradicated from the population [6]. When the bifurcation parameter ς = 180, then the mean time to loss immunity isT ≈ 90 days and the period of the oscillations increases monotonically with time, while the amplitude decreases with time and converges to an endemic equilibrium, see Figures 1(b ,d,f) and 2(b). When the bifurcation parameter ς = 2000, the mean time to loss immunity isT ≈ 1000 days and Figures 1(a,c,e) and 2(a) show waves of epidemic, with period and amplitude increasing with time until converging to a periodic solution of period T ≈ 500 days. It also appears that considering a given period, the cumulative number of infected individuals in the case ς = 180 is very high when compared with the cumulative number of infected individuals in the case ς = 2000 (figures not shown). This means that a vaccine or genetic selection of tilapia, aiming to increase the duration of waning will reduce the severity of the disease even if it will not completely eliminate the disease in the population. In both cases, the value of the basic reproduction number R 0 remains unchanged. This is explained by the fact that the expression of R 0 does not depend on the waning rate γ (·).
In contrast to the model developed in [19], the present model: • is able to reproduce all the scenarios developed in [19], namely the role of routes of transmissions of the TiLV. Therefore, it extends the model developed in [19]. • is able to study the influence of waning on the dynamics of transmissions of the TiLV, in particular, the waning plays a crucial role on the kinetics and the level of endemicity of disease propagation because the loss of immunity creates another time delay that can destabilize dynamics as shown in [2]: • when the temporary immunity period is within a certain range (about 1000 days), there will be periodic outbreaks of epidemic, and the disease will not be eradicated from the population, see on Figures 1(a,c,e) and 2(a) as shown in [6]. • when the temporary immunity period is about 90 days, there are some oscillations on the patterns, see on Figures 1(b,d,f) and 2(b). The existence of those oscillations was also obtained in [19], but their causes were unknown compared to the current model. • we observe that even if the basic reproductive number R 0 of both models are very closed (5.2630 and 5.2781), the disease dynamics in both cases are very different. • presents higher complexity and challenges to existing tools, some interesting results as the global stability of the disease-free equilibrium and the uniform persistence of the disease mathematically established in [19] were obtained here through numerical simulations.
The model developed in this paper relates to earlier work on modelling diseases with immunity. For example, Gomes et al. [16] have analysed several ODE models with temporary and partial immunity, as well as vaccination, where immunity wanes at a constant rate. They have shown that when a basic reproduction number exceeds unity, the solutions either decay linearly to an endemic steady state, or approach it in an oscillatory manner. In their findings, the immunity waning rate plays an important role in the time scale for oscillatory behaviour. While our model also highlights the importance of the waning and shows oscillatory behaviours of period depending on the waning rate. Since, temporary immunity is likely to wane at an increasing rate, we assumed that the duration of the immunity loss decays linearly with time since recovery. This led us to a simple form of the rate γ (τ ). However, with more biological data, one could adopt more complicated forms of the waning rate as in [5]. Indeed, in reality, the immune system of individuals may be boosted through exposition to the disease. This feature is a factor considered in several existing models through the immunity clock, that is by resetting the recovery age [27] or through the inclusion of additional internal states (within-host dynamics), see for instance [15,17].