key: cord-0007136-a4oten6u
authors: Korobeinikov, Andrei; Archibasov, Aleksei; Sobolev, Vladimir
title: Multi-scale problem in the model of RNA virus evolution
date: 2016-06-03
journal: J Phys Conf Ser
DOI: 10.1088/1742-6596/727/1/012007
sha: 89f8f963c728e580663a6eccdabaf49dcaa977ff
doc_id: 7136
cord_uid: a4oten6u

A mathematical or computational model in evolutionary biology should necessary combine several comparatively fast processes, which actually drive natural selection and evolution, with a very slow process of evolution. As a result, several very different time scales are simultaneously present in the model; this makes its analytical study an extremely difficult task. However, the significant difference of the time scales implies the existence of a possibility of the model order reduction through a process of time separation. In this paper we conduct the procedure of model order reduction for a reasonably simple model of RNA virus evolution reducing the original system of three integro-partial derivative equations to a single equation. Computations confirm that there is a good fit between the results for the original and reduced models.

Due to very high mutation and replication rates combined with high recombination abilities, RNA viruses are able to evolve very fast. This high evolvability makes RNA viruses a convenient experimental model in evolutionary biology and, at the same time, presents a serious problem, as viral evolution is probably the most important single factor accountable for the emergence of new pathogens and the development of drug resistance by existing ones.

Despite its apparent relevance, so far the development of mathematical or computational models of viral evolution which can be used in combination with experimental studies in evolutionary biology, remains a challenge even for such a simple object as RNA virus. In order to be useful in biological research, a mathematical model should be based on explicitly stated basic principles, postulates and hypotheses and allow a transparent interpretation of both results and parameters. In other words, a useful model should be essentially mechanistic. This implies that in evolutionary biology such a model must include a combination of factors responsible for natural selection (which drives evolution) and a mechanism for describing evolution. However, the factors generating the selection pressure act on a time scale comparable with an average generation time, or a life span of a single entity, whereas evolution is a comparatively slow process of accumulation of small changes that progresses through many generations to be evident. As a result, the model unavoidably includes multiple time scales which differ by several orders of magnitude. Computational or analytical handling of such models usually is a highly non-trivial task. In the case of viral evolution this discrepancy is further complicated by the fact that a typical generation times of virus and the hosts differ by orders of magnitude as well.

Besides phenomenological models, such as models suggested by Tsimring, Levin and Kessler [1] or Sasaki [2] and Sasaki and Haraguchi [3, 4] , are free of the mentioned drawbacks. However, the use of phenomenological models as a compliment to experimental research is usually questionable, as such models are rather illustrative than explanatory, and interpretation of models' parameters and obtained results in connection with an experimental layout is problematic.

These arguments imply that an ideal model should be mechanistic, include a number of properties of a real-life system and, at the same time, it should be sufficiently simple to allow analysis.

A relatively simple mathematical model of RNA virus evolution was recently suggested by Korobeinikov and Dempsey [5] . This mechanistic model is an extension of Nowak-May HIV model [6] , where viral phenotypes (strains) are assumed to be distributed in a continuous phenotype space, and random mutations are described by dispersion. For the sake of simplicity, Korobeinikov and Dempsey formulated the model on a basis of the 2-dimensional Wodarz-Christensen-Thomsen, or Wodarz, model [7] . The latter is a reduction of the original 3dimensional Nowak-May model under an assumption that the population of free virus particles is proportional to the infected cells population. However, an analytical study of even this comparatively simple model still is a highly challenging task, and it can be expected that for a more complex model, which includes factors, such as immune response or a therapy, the analysis would be even more difficult.

As we mentioned above, to a large extent the model complexity arises as a result of the presence of multiple processes progressing at very different time scales. Thus, the original Nowak-May model combines three processes and three time scales, namely (i) the proliferation of the uninfected target cells, (ii) the infectious process and infected cells' life cycle, and, finally, (iii) free virus life cycle. In the Wodarz model, which was used as a basis for model in [5] , one of these processes, namely the free virus life cycle, is omitted and replaced by the assumption that the free virus population is proportional to the infected cells population. The latter makes the model slightly simpler. However, model in [5] includes evolution which occurs on a comparatively slow timescale. Accordingly, the model is a slow-fast system. The considerable differences of the timescales, which make the model difficult for both computer simulation and analysis, at the same time, indicate that the model can be simplified using the process known as time scales separation (for instance, see [8, 9, 10] ). In this paper we carry out such a scale separation. As a result, the model RNA virus evolution in [5] is reduced to a single equation. Simulations demonstrate that apart from a comparatively short transition period there is a good fit of results for both, the full and the reduced models.

We would like to stress that the result is not limited by the specific model introduced in [5] . As the presence of multiple time scales which differ by a several orders of magnitude is typical for mechanistic models in evolutionary biology, the same technique can be applied to any such model.

As a basis for modelling, we employ the Nowak-May model of HIV dynamics [6] . This model describes interactions of three populations, namely uninfected target cells (for HIV, the target cells are T helper, or Th cells), infected target cells, and free virus particles, of concentrations u(t), V (t) and X(t), respectively. The model postulates that there is a continuous influx of the target cells (from the thymus, where they mature) at a rate b, that free virus particles infect target cells at a rate αXu, that the infected cells produce free virus particles at a per capita rate k, and that average life spans of the uninfected cells, infected cells and the free virus particles are 1/q, 1/m and 1/c, respectively. Under these assumptions the model equations are

The global qualitative behavior of this model is well studied and is known to be completely determined by the basic reproduction number of the virus (or that of infected cells, as for this model these numbers coincide) R 0 = αbk/cmq [11, 12, 13, 14] . Let us assume now that multiple viral strains are possible and are distributed in a continuous phenotype space. For model (1) a viral variant is characterized by parameters α, k, c and m, and hence it could be assumed that the phenotype space is 4-dimensional. However, it is known (and we will demonstrate this further) that for the Nowak-May model (1) these parameters are not independent: as we mentioned, for this model the viral dynamics is completely characterized by the basic reproduction number R 0 which usually serves as a measure of the Darwinian fitness of the virus. In other words, in the framework of the Nowak-May model a viral phenotype is uniquely characterized by the single number R 0 , and, therefore, it suffices to consider a 1dimensional phenotype space M = {s ∈ (0, ∞)}. Then fitness R 0 is a function of variable s. (The graph of this function is usually referred to as a "fitness landscape" . In the Nowak-May model coefficients α, m, k and c are attributed to a phenotype, and hence these are function of s. We assume that α = α(s) = α 0 α 0 (s), m = m(s) = m 0 m 0 (s), k = k(s) = k 0 k 0 (s) and c = c(s) = c 0 c 0 (s), where α 0 ,m 0 , k 0 and c 0 are coefficients of a wild strain (or a phenotype at the beginning of evolution), α 0 (s), m 0 (s), k 0 (s) and c 0 (s) are nondimensional functions and m 0 (s), k 0 (s), c 0 (s) = 1 for this wild phenotype. These assumptions lead to the following system of equations: and variable s is non-dimensional. Furthermore, u(t) is measured in cells · mm −3 ; v(t, s) and

x(t, s) are the densities of cells and free virus particles in the phenotype space, and hence these are measured in cells · mm −3 and virions · mm −3 ; cell production rate b is measured in cells · mm −3 · day −1 ; per capita mortality rates q, m, k and c are measured in day −1 , infectivity α is in mm 3 · virions −1 · day −1 , and dispersion coefficient µ is in day −1 .

A model of such type (based on the 2-dimensional Wodarz model, which is, in turn, a reduction of the 3-dimensional Nowak-May model) was introduced in [5] , and simulations with this model revealed a number of practically relevant outcomes. In particular, the model demonstrates the behavior which closely resembles the observed development of HIV infection. However, an apparent drawback of this model is that even in its present form its analytical study is a challenging task. A further development of the model, which could involve an incorporation of factors such as immune response and a therapy, would make the model more complicated and its analysis even more difficult. The complexity of this model is a result of an interaction of multiple processes with very different time scales. Thus, the original Nowak-May model (1) combines three processes and three time scales, namely (i) the proliferation of uninfected target cells, (ii) the process of infection and infected cells' life cycle, and, finally, (iii) the free virus particles life cycle. Compared with the Nowak-May model, model (2) includes the process of evolution, which, in comparison, is a very slow process, and hence model (2) is a slow-fast system. The presence of four considerably different timescales indicates that model (2) can be significantly simplified if we separate the timescales.

Let us denote t = Tt, s = Ss and u (Tt) = Uū(t), v(Tt, Ss) =Ṽv (t,s) , x(Tt, Ss) =Xx (t,s) , and assume that

Substituting these into system (2), we obtain a non-dimensional system

where ε = µ/q and ν = q/c 0 . Initial and boundary conditions for this system are, respectively,

and ∂v(t, 0) ∂s = 0, lim s→∞v (t,s) = 0.

Both numbers, ε = µ/q and ν = q/c 0 , are small parameters; for instance, for parameters, which were used in [5] and which correspond to a real life HIV case from [15] , ν ≈ 2 · 10 −3 and ε ranges from 5 · 10 −5 to 5 · 10 −3 . That is, system (3)-(5) is a singularly perturbed ("slow-fast") system, and its size can be reduced using the scale separation technique [8] . This procedure corresponds to assumptions that (i) there is a stable slow manifold such that all solutions converge to this manifold, and hence (ii) the values of fast variables u(t) and x(t, s) converge to quasi-equilibrium values which slowly vary as the system progresses along the slow manifold [9] . Let ν = 0, then equalityk (s)v =c(s)x holds. Substituting this equality into (3) and (4), we obtain equations

where β(s) = α(s)k(s)/c(s). System (6), (7) was suggested in [5] . Now let ε = 0, thenū

whereλ(t) = (7), we finally obtain a single equation

where R 0 (s) = bα(s)k(s) qm(s)c(s) ,m(s) = m(s) µ .

When λ(t) q, that is when the reduction of the uninfected population due to infection is small, equation (8) can be further simplified: in this case we can employ the Taylor expansion and use the approximate equality

Then 

Following [5] , and for the sake of simplicity, we assume that the diversity of viral phenotypes is described by a single parameter α(s) which can be interpreted as the efficiency of a single virus particle in infecting a target cell; the other parameters are assumed to be constant and the same for all genotypes. To some extent this assumption is justified by the fact that for the Nowak-May model genotypes fitness is defined by a single number R 0 = αbk cmq . Under this hypothesis, and taking into consideration that R 0 is a function of all six parameters, it should make no difference which of these six parameters is varying in the phenotype space. The choice of α(s) is just a matter of convenience, because the direct dependence of R 0 on α simplifies the notation.

Results of numerical analysis are presented at the following figures. (For convenience of comparison, results in these figures are shown for physical dimensional variables.) For these simulations we use the following values for the system parameters: b = 20 cells · mm −3 · day −1 , m = 0.8 day −1 , q = 0.02 day −1 , k = 10 4 day −1 and c = 8 day −1 . These parameters correspond to patient 2 in [15] . For simplicity we assume a linear landscape, postulating that α(s) = as.

Here a = 8 · 10 −7 mm 3 · cells −1 · day −1 ; this implies that R 0 = 1.25s. All simulations are for µ = 10 −6 day −1 .

Simulations demonstrate an excellent fit of results obtained for system (6), (7) and equation (8) . shows that results for these two systems perfectly coincide everywhere apart from a comparatively short transition period. (We have to remind that equation (8) is unable to describe the transition dynamics at all.) If the initial conditions are taken on or near the slow manifold, then the results for these two systems coincide everywhere (see Fig. 2 ). However, a match of results for systems (3)- (5) and (6), (7) (and hence for equation (8)) is only qualitative. Figure 3 demonstrates that system (3)-(5) exhibits a slower evolution, producing outcomes which are delayed in time, compared with these for system (6), (7), or equation (8). We suppose that this effect is caused by a continuous accumulation of a small delay, caused by the third equation of the Nowak-May model (2) , over a long period of evolution.

Due to the fact that a mathematical model in evolutionary biology should necessary combine the factors, which characterize natural selection and drive evolution, and which acts on a comparatively fast time scale, and a very slow process of evolution itself, such a model is complex and its analysis is difficult. However, a combination of several time scales within a model suggests that such a model can be significantly simplified using the scale separation techniques. For ODE such technique is based on the famous Tikhonov's theorem. The mathematical models in evolutionary biology are usually formulated as integro-differential equations and PDE. However, with some reservations the same concepts can be applied to these equations as well (for instance, see [10, 16, 17] ). In this paper we execute the scale separation for a RNA virus evolution model which is based on the 3-dimensional Nowak-May model. As a result of this procedure, the original system of three integro-partial differential equations is reduced to a single integro-partial differential equation, describing the dynamics of infected target cells v(t, s). The other variables, namely the uninfected cells u(t) and free virus particles x(t, s), are expressed as explicit functions of variable v(t, s). Numerical simulations show an acceptable fit of the solution of the original and the reduced models.

A specific feature of this particular model, which is based on the Nowak-May model specifically tailored to describe the HIV dynamics, is that target cells are assumed to be immigrating from the thymus at a constant rate. However, the same concept and the same techniques can be applied to a model of evolution based on any other model of virus dynamics.

RNA virus evolution via a fitness-space model

Evolution of antigenic drift/switching: continuously evading pathogens

Evolutionary pattern of intra-host pathogen antigenic drift: effect of crossreactivity in immune response

Antigenic drift of viruses within a host: a finite site model with demographic stochasticity

A continuous phenotype space model of RNA virus evolution within a host

Virus dynamics: Mathematical principles of Immunology and Virology

The importance of lytic and nonlytic immune responses in viral infections

Stability analysis of nonlinear multiparameter singularly perturbed systems IEEE

Singular Perturbation and Hysteresis

it The Boundary Function Method for Singular Perturbation Problems

Global properties of basic virus dynamics models

Global asymptotic properties of virus dynamics models with dose dependent parasite reproduction and virulence, and nonlinear incidence rate

Stability of ecosystem: Global properties of a general prey-predator model

Global stability of a population dynamics model with inhibition and negative feedback

Modeling plasma virus concentration during primary HIV infection

Global region of attraction of a periodic solution to a singularly perturbed parabolic problem

The Cauchy problem for a singularly perturbed integro-differential Fredholm equation