From ODE to Open Markov Chains, via SDE: an application to models for infections in individuals and populations

Abstract:We present a methodology to connect an ordinary di erential equation (ODE) model of interacting entities at the individual level, to an open Markov chain (OMC) model of a population of such individuals, via a stochastic di erential equation (SDE) intermediate model. The ODEmodel here presented is formulated as a dynamic change between two regimes; one regime is of mean reverting type and the other is of inverse logistic type. For the general purpose of de ning an OMCmodel for a population of individuals, we associate an Ito processes, in the form of SDE to ODE system of equations, by means of the addition of Gaussian noise terms which may be thought to model non essential characteristics of the phenomena with small and undifferentiated in uences. The next step consists on discretizing the SDE and using the discretized trajectories computed by simulation to de ne transitions of a nite valued Markov chain; for that, the state space of the Ito processes is partitioned according to some rule. For the example proposed for illustration, the state space of the ODE system referred – corresponding to a model of a viral infection – is partitioned into six infection classes determined by some of the critical points of the ODE system; we detail the evolution of some infected population in these infection classes.


Introduction
In this work we have two generic goals. The rst goal is to show how to associate an initial ODE model applied to a generic individual of a population to a Markov chain model for the evolution of an open population of these individuals. The second goal is to illustrate the details of the procedure by considering the case of a viral infection on a human population.We thus show how to build a nite state Markov chain model starting from an ordinary system of non linear di erential equations. The method relies on an intermediary step that associates a system of stochastic di erential equations to the ODE system.
In [15] there is a reference to the possibility of associating a Markov chain to a stochastic di erential equation with the ultimate purpose of studying a model, by simulation, in a similar manner as the well known Markov Chain Monte Carlo method. The study of the properties of a stochastic di erential equation of di usion type that transfer to an associated Markov chain, by means of a discretization scheme -such as the Euler-Maruyama -has been made almost two decades ago, see [33]. Other important references that study properties of an associated Markov chain, created from a di usion by means of a discretization scheme, are [31], [29] and [19].
There are several examples in the literature of the type of models studied in this work, namely, identifying particular regimes. For instance, in [4] the authors focused on the mean reverting type of model for the SDE of interest. In [35], using Kalman lter approach for the parameters estimation the authors consider a di usion with constant volatility and linear drift.
In [34] coupled systems of two and three SDE with constant volatility are considered. In [9] a modi ed Gompertzian SDE equation, with linear volatility function and an explicit solution similar to the Ornstein-Ulhenbeck process, is used to describe bacterial growth under bactericide. In [25] the approach is also from an ODE model to a SDE system model in order to account for the fact that (and we quote) process is not smooth, subject as it is to a variety of metabolic and hormonal in uences, which change over time. The system has only one equation, with volatility term proportional to the product of the two processes and is coupled as the two processes appear in the equation with null volatility.
The most used estimation methods for SDE do not apply to the model here proposed, see [5] for a comprehensive review on the subject. We anticipate several di culties developing such methods and we hope to address this issue in future works. For the moment we parametrize the SDE model to obtain illustrative examples for which the methodology here proposed can be applied.
The importance of linking within-host dynamics and between-host dynamics has been long recognized, see [14]. Heterogeneous population dynamics results from individual di erences which produce the complex epidemiological patterns observed. For virus infections, in particular viral infection, authors have considered e.g. the connection between the two scales through the epidemiological reproduction number [3] and, using two ODE nested models to study the evolution of virulence [10].
Let us comment on some aspects of the example of viral infection we chose to illustrate our method. It is well known that there is an obvious relation between the number of viruses attacking leucocytes. As the viral load increases the leucocyte count is expected to decrease; moreover we should expect a regular variation in this coupled evolution. Being so, an ODE model for an average behavior seems appropriate. But, as there are many other factors conditioning the evolution of the mutual in uence process, it is also expected that, for each individual data set variations on the average evolution pattern may appear. In [26] there is a reference to in uences, in a particular process interfering with the interaction of viruses and leucocytes, that are classi ed as being stochastic in nature. Being so, a SDE version of the initial ODE model for the pure virusesleucocyte interaction may be justi ed. Some other work on the same lines of thought is the analysis of a three dimensional model of virus competition that is presented in [32]; a ve dimensional analysis, contemplating both infected and non infected leucocyte cells is given in [6]. Let us brie y describe the contents of this work.
• An overview of the methodology proposed in this work is presented in Section 2. • In Section 3, the very important idea of competing regimes is presented and next we model the coupled evolution of two competing entities in an individual -for instance, a viral load and a leucocyte countby a system of two ordinary di erential equations. This ODE model is supposed to describe the average behavior of the quantities modeled. • Secondly, in Section 4, we consider a stochastic di erential equation model (SDE) built upon the ODE model by adding a noise term to each of the ODE. Each trajectory of this SDE models the evolution of the two competing and interacting entities for a particular individual. As usual, the interpretation of a SDE model built by adding a noise term to the equations of a ODE model is that we consider that the evolution of the modeled quantities for an individual -which in general is described by the ODE model -is perturbed, at the individual level, by unknown factors. • The third step -developed in Section 5 -consists on associating a Markov chain transition matrix to the SDE model, by simulating trajectories of the solution process of the the SDE and counting the transitions between Infection Classes. These Infection Classes were de ned by regions in the phase space of the ODE by means of critical points of the solutions of the ODE system. • In the last step -presented in Section 6 -we consider an open Markov chain (OMC) model for the population with ve transient states (the rst ve stages of infection) and an absorbing state corresponding to the last stage of infection, from where we assume that it is not possible to move out without treatment or death. The OMC model is adequate since the population is open and constantly fed by the number of newly viral infected. For general conditions on the number of newly viral infected patients, we are able to determine the evolution of the expected numbers of patients of the subpopulations -the Poisson parameters -in each one of the Infection Classes -as done in [7]. We exploit the open Markov chain model by considering simulated data from a yearly viral newly infected population on a 23 year period.
Remark 1. As so, this sequence of models, each one built upon the preceding one -ODE, SDE and OMC models -allow to transfer information, in a coherent way, from an individual level to a whole population level.

From ODE to Markov chains via SDE: an overview
The general pathway from an ordinary di erential equations model to open population Markov chain model via a stochastic di erential equation model is schematically described in the following. Let us suppose that we have a model given by the solution of a Cauchy-Lipschitz ordinary di erential equation (ODE) -satisfying some regularity assumptions -and an initial condition such as: or, integrated form, Now, suppose that we consider equation (2.1) as giving the evolution of the mean value of a population size and that, aiming at having an individual random model, we add a noise, given by an Ito integral with respect to a Brownian process (B t ) t∈[ ,T] -de ned on a probability space (Ω, F, P) -such that we get a stochastic di erential equation (SDE) given, for ω ∈ Ω, by: f (s, y(s, ω))ds + t t σ(t, y(t, ω))dB t (ω) , (2.2) or in the usual di erential form, d y(t) = f (t, y(t))dt + σ(t, y(t))dB t , y(t ) = y .
Let us observe that under su cient hypothesis guaranteeing that the Ito stochastic integral part of equation (2.2) is a martingale, with constant null mean, we will have, From this point on we will not distinguish y from y. Now, for illustration purposes only ¹, suppose that we discretize equation (2.2) with the Euler-Maruyama discretization scheme to get: y n+ = y n + f (tn , y n ) (tn+ − tn) + σ(tn , y n ) (Bt n+ − B tn ) , for t < t < · · · < t N = T a subdivision of [ , T]. In the unidimensional case it is known that, and we have that the sequence (Un) n= ,...,N− , de ned by: with F ← n the inverse of the distribution function of a random variable distributed with the law N ( , (tn+ − tn)), is a sequence of independent uniformly distributed random variables. In this way, if , that is, the di usion given by (2.2) is homogeneous, with, we may write that, y n+ = G(y n , Un) , n = , . . . , N − , which is the functional de nition of a Markov chain (see [28, p. 62]). As so, we could try to apply the methodologies used in [7] and [8] -and other references therein -with the purpose of identifying the possibility of stable behavior in the set of transient classes. Let us observe that, in principle, we will have that the chain (y n ) n= ,...,N , for large N, may have a large number of states and so, in order to apply the referred methods of open Markov chain models, a natural situation would be to de ne a small number of classes C i by considering a partition of the state space, S i , with S = ∪ d i= S i and S i ∩ S j for i ≠ j, such that at step n, ω ∈ C i ⇔ y n (ω) ∈ S i .
In the context of epidemiological studies this is feasible if the state space is some set of values of some biological indicators or markers like a virus load and some cell count.

A non linear ODE 2-dimensional model
The competing regimes modeling we adopt in this work is rooted in homeostasis which is the maintenance of an internal stable environment for proper functioning of the body. As a consequence, we can suppose that body markers are subjected to variations that tend to conform to some alternating regimes, the alternation being driven by other processes in the body. For simplicity we suppose that, for the body marker we are going to model, the values can evolve in two di erent regimes. For illustration purposes we consider the usual mean reverting model given by the ODE, which supposes that the marker has a tendency to evolve towards the value 1; and the inverted logistics model given by, which essentially says that the marker will evolve towards some larger value. We present some solutions of these ODE, corresponding to di erent initial values in Figure 3.1. If in the homeostasis context one kind of evolution is changed in a di erent kind of evolution by the e ect of external causes, we may consider that these two regimes are coupled by a regular function that shifts, in time, the evolution from one regime to another; the coupling will be provided by a convex combination of the two right-hand members of equations (3.1) and (3.2). For technical reasons we suppose that the coupling function, depends on two parameters α and λ ² and is given by: Solutions of The ODE with the two regimes, coupled by F, is given by (3.4) with additional parameters µ controlling the intensity of the reversed logistics drive, uMin, which is a bound on the values of u and ν, which is the mean reverting level.
In Figure 3.2 we show, both F, for λ = and for several values of the parameter α and a solution u(t), t ∈ [ , ] -for α = . in the coupling function f -of the ODE with the coupled regimes with some additional parameters given in Table 1.   For modeling purposes we may use this particular kind of evolution, the one described by u(t), t ∈ [ , ] in the right-hand side of Figure 3.2, for some biological internal marker of a normalized count that evolves with a continuous transition between a decreasing regime and an increasing one, according to the variation of the function f . We observe that we could also have a continuous transition between an increasing regime and a decreasing one by changing the weights f and − f , in the right-hand side of the ODE.
Remark 3. We stress that further combinations of regimes are possible: either by considering other coupling functions performing a convex combination of more than two terms or by considering a second convex combination -with another coupling function -of the existing convex combination. This could be done to alter the asymptotic regime. Also, we observe that the parameter α de nes, in some sort, the speed of the regime change.
There is a natural perspective to study epidemics namely, the population dynamics angle (see, for instance, [18] and references therein). To study the infection at the individual level and on the average, the approach relies on the modeling of the dynamics -by means of ODE -of some body internal markers, such as viral charges or leucocyte counts (see [1]) or leucocyte percentages or ratios.
That is the general approach we pursue in the following considering the rst part of the proposed methodology: a non-linear ODE model for infections in which we may have two interacting entities -for instance some viral load and for instance the leucocyte count. This model, which aims at describing an average general behavior, has two coupled ODE, one for u -the leucocyte count -and another for zviral load. In each equation there are two evolution regimes and the transition in time, between these regimes, is achieved by a coupling function. We stress that the model chosen is not of mechanistic type in the sense that it does not explicitly model the interaction of the two populations as in [23] and [27].
Let us detail some general features of this model: 1. Both the equations have two regimes, of the same type in the two equationsinverse logistics and mean reverting -coupled by the coupling function F(t, α, λ); this function of logistic type goes from zero -at time zero -to values very close to one in nite time. 2. The initial regime for u, controlled by − F(t, α , λ ), for small values of t, is the inverse logistics term for the function z, which is also the nal regime, controlled by the values of F(t, α , λ ) for the function z.
Being so, in this model, the viral load exerts some in uence in the initial behavior of leucocyte count.

. The qualitative study of the ODE model
In this section we state some questions related to the qualitative study of the set of ODE given in (3.5). We observe rst that this system can be seen as an autonomous system by considering: (3.6) In fact, let us consider, for α ∈] , [, the di erential equation: and so the equation in (3.7) may be written as with the immediate solution v = − γe −λ( −α)t , with γ > some integration constant, which gives, as the general solution of equation (3.7). The qualitative study of a system such as the one in (3.6) -a system of four non linear ODE -may require some special delicate analysis. An alternative path could be to study directly the system (3.5) by means of some ad-hoc procedures. Examples of such studies are [11], [17] and [2].
Remark 4. The linearized part of the ODE system (3.6) is given by: The linearized part of the ODE system (3.6) when u = z = r = s = is given by: The eigenvalues of the linearized part at ( , ) of the system given in (3.6), given as function of the parameters, are the following.
The rst and the third ones are negative and the second and the fourth ones are positive. As so, by the theorem of Hartman-Grobman, see [24, pp. 120-127], the trajectories of the linearized system are homeomorphic images of the original system in a neighborhood of and so, the point ( , ) is not a stable point of system (3.6). This property may have relevance for the question treated in subsection 10.

. An illustrative numerical application
The system of equations (3.5) was solved numerically using the parameter values given in Table 1. These parameters were chosen in order to provide a solution with the desired characteristics.
Remark 7. For practical application of this model, the estimation of the parameters for equation system (3.5) should be done using the corresponding set of Ito processes in unidimensional case, in order the methodology to be applied with a set of signi cant data. E cient estimating methods are being developed. Also, the qualitative study of the system of equations (3.5) should provide an insight into the desired properties of the model and the corresponding range of the parameters.
In Figure 3.3 there is a representation of the solution (u, z) of the system of equations (3.5) with parameters given in Table 1. This solution was obtained by numerical integration ³.   Figure 3.3. 2. The present report may be thought as a proof of concept exposition so, no homogeneity analysis for the units choice has been performed. We must stress that, in principle, the values of the viral load must be thought to be given in log scale. At this point of the exposition the time scale presented was not adjusted to any real data; see Remark 19 for an adjusted time scale proposal. 3. The pattern of the initial decreasing of the viral load that appears in Figure 3.3 was rst modeled in the works of M. Novak and R. May. See, for instance, [22] where several models, with complex dynamics explaining this pattern, are studied.
We must emphasize that the model here studied does not take into consideration the existence of some kind of disease treatment; for such a purpose, within the approach developed in this work, we would add a third regime describing the average response of both variables to such a treatment. One example, within a di erent context, of such model is presented in [16].

The associated SDE 2-dimensional model
We now develop the second part of the proposed methodology: to model the two entities interaction -for instance, viral load and leucocytes count -at the individual level, by adding a Wiener process type noise to the ODE model introduced in Section 3. The SDE system we consider is given by: where (B t ) t≥ and (B t ) t≥ are two independent Wiener processes and where, for notational simplicity for the coupled system of SDE, we have de ned the drift functions taken from the right-hand side of the ODE system in formulas (3.5), by: , Remark 8. We note that the additive noises used -having an intensity proportional to the square root of the SDE equations unknowns -X t (respectively Y t ) , corresponding to u (respectively v) in the ODE model -were chosen in order to ensure, if possible, that the solutions of the SDE system are always non negative.
Remark 9. Although it would be possible to consider some non identity xed covariance matrix between (B t ) t≥ and (B t ) t≥ , in formulas (4.1), we have assumed independence by lack of an interpretation of this possible covariance structure. This covariance structure may be added, in a real data model, upon the development of the parameter estimation procedure.
The system of SDE given by formulas (4.1) was discretized in order to obtain a simulated sample of trajectories.
The parameter values used for the simulation are shown in Table 2. The trajectories of both processes de ned by the SDE system (4.1), (X t ) t≥ and (Y t ) t≥ , which were obtained by the simulation procedure, can be seen in the left-hand side of Figure 4.1. Remark 10. A far easier approach to the conjectured result on positivity of the kind of di usions used in this work is to show that the discretized sequence will only take positive values. This amounts to use properties of Gaussian random variables and can be done by induction. The stability analysis of the ODE system (3.6), in Subsection 3.1, has consequences for the proof of a technical result ensuring that the solutions of the SDE system (4.1) will remain positive almost surely. The general idea for such a result is that as the processes approach zero, with positive values, the volatility part of the equation becomes less and less important and so the drift part of the equations -that is the part corresponding to the ODE system given by (3.5)-dominates; as the point ( , ) is not stable (see remark 6), by force of two positive eigenvalues of the linearized part of the SDE system, the solution should be forced to stay on the positive values. The papers [21] and [20] contain recent work that may be a starting point to the necessary technical proof on the positivity of trajectories of the SDE of the model here presented.  Simulated trajectories for X t , the viral load (orange), and for Y t , the leucocyte count (blue), SDE model and, on the right, the parametric plot of (u(t), z(t)) together with the images by u (respectively z) of the critical points of u, u (respectively z, z )

The associated Markov chain model
We now detail the procedure for the third part of the proposed methodology: associating an OMC model to our SDE model. We follow this procedure in two steps: Step 1 We partition the state space of the SDE in six regions, corresponding to six Infection Classes, by means of the critical points of the solutions u and v of the ODE system (3.5) (see Figure 4.1 for a quick identi cation of these regions and Table 4 for the precise de nition of these regions). The determination of a set of critical points of u, u and z, z and the corresponding images by the functions u and z, are presented in Table 3. These points correspond, in the right-hand side of Figure 4.1, to a vertical line -for u(ta) -  and two horizontal lines -corresponding to z(t b ) and z(tc) -allowing for the de nition of six regions in the phase space of (u, z) that are indicated in Table 4. These regions admit qualitative descriptions as

z(tc)[ [u(ta), +∞[×[z(tc), +∞[ [ , u(ta)[×[z(tc), +∞[
Infection Classes. Region 6 is the terminal infection state with a higher viral load and a reduced leucocyte count. Region 4 corresponds to the infection establishment followed by a noticeable increase of the viral load and the corresponding immune response. Region 2 is the most reasonable initial infection state due to its de nition in the phase space of the SDE and so on and so forth.
Step 2 We simulate trajectories of the SDE (4.1) and we count the transitions between the di erent regions corresponding to the six Infection Classes. That is, considering every pair of regions and by counting the transitions from the rst region of the pair to the second region of the pair, for every trajectory of a sample -of (only)16 -simulated trajectories, we get the transition matrix T presented in (5.1). Increasing the number of trajectories does not noticeably change the accuracy of the estimation of the Markov chain transition probabilities.
We stress some general methodological remarks for the building up of the Markov chain model associated to the SDE extension of the ODE model: 1. The determination of the state space partition for the determination of the states of the Markov chain model could be done, a priori, by using, rstly, the critical points of u and v and after the in ection points, that is the critical points of u and v . The interpretation of (u, v) as the mean value functions of the SDE solution trajectories may give a justi cation for this choice. The computation of these critical points is feasible after the parameters are estimated. This is the approach in the following developments.  3. The simulated trajectories must be determined in a way that preserves the -possible -ergodicity, and other asymptotic properties -of the di usions discretized solutions of the SDE model.
A quick analysis over matrix T, reveals that the Markov chain corresponding to this transition matrix is irreducible and therefore all states are recurrent. So, it may be that the initial proposition of considering the partition of the state space using the critical points is not adequate in order to obtain a matrix with transient states and possibly one only recurrent state. In any case, the partition of the state space should be performed according to some justi cation drawn from the phenomena under study.
Remark 11. Besides the classes de nition, there are other choices that may have an in uence on the values of the transition matrix; these are the discretization step -which in this case was taken as 0.01 -and the discretization scheme chosen. As already said, the choice of the discretization scheme for the SDE has an in uence on the properties of the di usion that are transmitted to the Markov chain, such as the ergodicity.

. The altered Markov chain model associated to the SDE system model
In this subsection, we will justify the assumption that the sixth state -corresponding to the sixth Infection Class -is the sole recurrent state and that all the other states are transient. This assumption, besides being justi ed, is also necessary for the OMC model implementation (see [12] for a full development of these properties). The assumption leads us to change the transition matrix T into T, in (5.2), by simply considering non existing the transitions from class 6 to classes 5, 4 and 3; we stress that in T these transitions have small probabilities, respectively .
Remark 12. We again recall that we are not taking into consideration the existence of disease treatments. However, even if we did so, there would always be a recurrent state, namely death, and so this highlights that an OMC model would be adequate to describe viral infected populations.
Remark 13. We could, justi ably, alter matrix T by considering that all entries smaller than a certain threshold -say . , for instance -would be considered zero and the corresponding probability would be summed either to the greatest probability in the row or distributed proportionally through the non zero entries of the row. Any of these choices should be done according to the interpretation of the classes and respective transitions. An experience performed during our study did not show any appreciable changes in the conclusions.
With this altered matrix we now have a Markov chain model with the rst 5 states being transient states and the sixth being a recurrent state. The evolution of the proportions of the elements of some closed population -an unique cohort -with initial distribution given by: { . , . , . , . , ., .} , is presented left-hand side of Figure 6.1 (see [8, p. 280] for a quick mathematical justi cation). As expected, the proportions in the transient classes -Infection Classes 1 to 5 -decrease while the proportion in the Infection Class 6 goes to one. This is a natural asymptotic behavior of an usual closed population Markov chain model.
The initial distribution (5.3) was chosen according to the following idea: in principle, in a country with a wide reaching health care system, the newly viral infected will be mostly detected whenever belonging to Infection Class 2; the one class where there is already a shortening of leucocytes and the presence of viral load. The results of the evolution, that we show, depend, of course on the initial distribution.
Remark 15. With the purpose of showing the quantitative e ects of the initial distribution, on the evolution patterns of the population, we considered a second initial distribution given by: { . , . , . , . , . , .} . (5.4) This corresponds to more weight in the Infection Classes 1, 3, 4 and 5 and corresponding less weight in Infection Class 2. The evolution dependence on the initial distribution highlights the need for an accurate estimation of this initial distribution.

Open Markov chain modeling: a review
We now detail and comment the results that will be used in this paper to get the evolution of the newly viral infected in the di erent Infection Classes.
The study of open Markov chain models we will present next relies on results and notations that were presented recently in [7] and that we reproduce next, for the readers convenience. We will suppose that, in general, the transition matrix of the Markov chain model may be written in the following form.
where K is a k × k transition matrix between transient states, U a k × (r − k) matrix of transitions between the transient and the recurrent states and V a (r − k) × (r − k) matrix of transitions between the recurrent states. A straightforward computation then shows that, We write the vector of the initial classi cation, for a time period i, as with t i the vector of the initial allocation probabilities for the transient states and r i the vector of the initial allocation probabilities for the recurrent states. With (6.1) and (6.2), we now notice that the vector of the Poisson parameters, for the population sizes in each state at time T, may be written as We observe that the rst block corresponds to the transient states and the second block, the one in the righthand side, corresponds to the recurrent states. From now on, as a rst restricting hypothesis, we will also suppose that the transition matrix of the transient states, K, is diagonalizable and so, with (η j ) j∈{ ,...,k} the eigenvalues, (α j ) j∈{ ,...,k} the left eigenvectors and (β j ) j∈{ ,...,k} the right eigenvectors of matrix K. We observe that j ∈ { , . . . , k} corresponds to a transient state if and only if | η j |< . We may write, and so, as a consequence of (6.3), for the vector of the Poisson parameters corresponding only to the transient states, λ +t T , we have: The main result we will use next is the following.
Remark 16. We observe that proportions in the Markov chain transient classes, on both statements of the Theorem 6.1, only depend on the eigenvalues η j , j = , . . . , k. In fact, whenever using formula (6.6) to compute proportions these proportions do not depend on the value of λ as we have that, and the term in the right-hand side multiplying λ is a vector with the dimension equal to the number of transient classes k, which is equal to the dimension of the square matrix K. As so, when computing proportions, by normalizing this vector with the sum of its components, λ ≠ disappears.

. . The main theorem formula interpretation
This subsection aims at clarifying some aspects of the implementation of the main theorem 6. In [30, p. 111] it is shown that if x , . . . , xm is any set of orthonormal eigenvectors corresponding, respectively, to distinct eigenvalues η , . . . , ηm of a matrix A then, Also, if the eigenvectors are not an orthonormal set, see [30, pp. 169, 170], if the matrix A is diagonalizable with eigenvalues η , . . . , ηm and the correspondent independent eigenvectors x , . . . , xm, then:

. The OMC model associated to the SDE system
In an OMC model we consider a permanent in ow of new members in the population which has an evolution driven by a Markov transition matrix. In an OMC model the asymptotic behavior of subpopulations in the transient classes become non trivial and may be quantitatively described. So, in order to study the evolution of the whole population of newly viral infected individuals spread in a long period of time -by an OMC model -we will consider, as an example, simulated data of newly detected with some infection on a period of, say, 33 years, for instance from 1983 to 2016. In the right-hand side of Figure 6.1 we present the tting of this data with a function such as, Remark 17. The parameters a, b, c, d, k for the best tting of the function g a,b,c,d,k given in formula (6.10), may be estimated -in a more rigorous way -by maximum likelihood, as done in [7] for a sigmoid type function or in [13] for an exponential type function.
Remark 18. The function g tted used to t the evolution of the newly viral diagnosed has a stationary point xstationary ≈ and after it becomes increasing. A natural alternative choice in the absence of further information and for illustration purposes lead us to consider the function G, which coincides with the tted function -see (6.10) -until this stationary point and remains constant for values of the variable greater than this stationary point.  As an application of Theorem 6.1, namely of formula (6.6) we obtain the estimates of Table 5 and Table 6. * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *   When the number on newly diagnosed increases without upper bound from date 70 -but under the condition stated in the second thesis of Theorem 6.1 -formula (6.7) shows that nevertheless the mean value of number of elements in each Infection Class, divided by mean value of the total newly diagnosed at each date, also stabilizes as a constant value. This can be seen in Table 6 where we present the correspondent results.  Tables 5 and 6. 1. As a consequence of the assumption made that the newly diagnosed viral infected come from a Poisson distributed independent populations in each year, the populations in the di erent Infection Classes 1 to 5, are also Poisson distributed and the evolution of the Poisson parameters for these populations -and so also the expected number of these populations -is what is shown in the right-hand side of Figure 6.2. 2. The rst obvious conclusion is that the evolution of the proportions in the transient classes corresponding to the di erent cohorts -dates 1983-2016 -is quite di erent for the proportions evolution of one only cohort as shown in Figure 6.1. 3. The evolution of population in Infection Class 2, clearly, re ects the best the evolution pattern of the newly diagnosed viral infected as shown in Figure 6.1. This a direct consequence of the initial distribution in which we assumed that % of the newly infected would be in this Infection Class 2. 4. As a consequence of Remark 16, the proportions of the elements in each of the Infection Classes are the same either when the newly diagnosed are in a constant expected number from date 70 on or when the expected number of the newly diagnosed grows without bound from date 70 on. And so, the bottom half of both Tables 5 and 6 are equal, despite being computed by averaging di erent vectors.