key: cord-1041367-nfhfys7r
authors: Bartoszek, Krzysztof; Bartoszek, Wojciech; Krzemiński, Michał
title: Simple SIR models with Markovian control
date: 2021-02-16
journal: Jpn J Stat Data Sci
DOI: 10.1007/s42081-021-00107-1
sha: f1788e594741d28300eb0ac96a10be0e0f8ac8c6
doc_id: 1041367
cord_uid: nfhfys7r

We consider a random dynamical system, where the deterministic dynamics are driven by a finite-state space Markov chain. We provide a comprehensive introduction to the required mathematical apparatus and then turn to a special focus on the susceptible-infected-recovered epidemiological model with random steering. Through simulations we visualize the behaviour of the system and the effect of the high-frequency limit of the driving Markov chain. We formulate some questions and conjectures of a purely theoretical nature.

The spread of an epidemic and its related characteristics in a large population may be efficiently described by deterministic models. In particular, differential (or difference) equations are a popular first choice. However, if the population (or a particular subpopulation which plays a specific role in the epidemic) is small, then a stochastic approach is more appropriate. For example, at the initial stage of an epidemic when the disease is brought in from outside the considered population, only single individuals may be infected or we may even have the extreme case of a single individual

(the so-called patient zero). Moreover, these infected individuals can be hidden from the health authorities. Our approach and results are general and can be also considered in epidemics amongst livestock (e.g. the recent spread of the African swine fever virus in Eastern Europe). We also do not focus on any particular mechanism of infection by the driving process. It can be due to transmission from natural interactions or even intentional contamination (from e.g. market competition). Alternatively, we may be in a situation, where the disease spread is impossible to observe. This is, in particular, true with the current dynamics of the SARS-CoV-2 virus. What is observed is the number of confirmed cases (based on tests) and confirmed case fatalities. However, there is a very large number of asymptomatic cases that are not discovered (Bai et al. 2020) . These people are contagious and this is where the dynamics of the virus are taking place. However, we can only model, e.g. via SIR (susceptible-infected-recovered) systems, those parts of the population for which we have data. Therefore, it might be more appropriate to consider an IR, random dynamical system (RDS) model (for an overview of RDS with multiple examples see e.g. Smith and Thieme 2011; Ye et al. 2016; Ye and Qian 2019) . The Infected dynamics are controlled by a random Markov process-i.e. from a randomly behaving population, we obtain observed infected individuals. There are two arguments for considering the unobserved dynamics to be random. Firstly, despite the population being large, various official or unofficial (e.g. self-isolation) mitigation strategies could result in random effects. This can be due to the randomness (from the perspective of the observer) in the moment of implementation, their scale, and also to what extent they are carried out by the population. Secondly, as one only tests a certain (small) part of the population and not completely at random but according to some protocol (e.g. only symptomatic, people who had contact with confirmed cases, risk groups, medical personnel, etc.), then one can expect randomness in the response-infected people randomly entered the particularly tested group. This is even true for the testing of symptomatic individuals-the symptoms of COVID-19 are consistent with many other illnesses, like the common cold or the flu. In all these situations the ignition points of the infection are limited in numbers. Therefore, one should not ignore all mathematical conditions of differential (smooth) models, like SI, SIR, etc. In this work we try to merge deterministic and stochastic methods in order to obtain a combined, hybrid model. We follow the approach of dividing the whole population into subgroups, arriving at a typical compartmental model. We assume that the total size of considered population is fixed and does not change through time. After normalization, we assume it to be simply 1. We follow standard nomenclature and call the healthy members who are exposed to the disease as susceptible (the total amount of the exposed is S), those who are infected are denoted by I, and finally, the "absorbing" group, the recovered (or removed from the population) is denoted by R. We have that S + I + R = N , where (in general) N ≡ const , does not depend on time, and as already mentioned we take N = 1.

The main aim of our paper is to illustrate the dependence of the quantitative analysis of an epidemic on a random steering process (which could represent the intensity of social relations). First in Sect. 2 we introduce the necessary mathematical background and our main results. Then in Sect. 3 we present two types of simulations. One, a toy model of a linear random dynamical system, and then a simple application, based on the SIR model to epidemiological data. In Sects. 4 and 5 we give a careful treatment of SIS and SIR models with our introduced random control mechanism We end our work with some some historical and outlook remarks, Sect. 6.

The SIR model is in fact a flow (S(t),

, governed by the system of differential equations. We wish to consider a wider setting by introducing general dynamical systems (cf. Smith and Thieme (2011) ).

is a semigroup, i.e. a representation of the semigroup ℝ + as transformations of X.

Definition 2 A semiflow is called state continuous if for every t ∈ [0, ∞) the mapping t ∶ X ↦ X is continuous. is called uniformly continuous if for every > 0 there exists > 0 , such that if (x, y) < , then is called uniformly continuous in finite horizon time if for every T > 0 and all > 0 there exists a T > 0 , such that if (x, y) < T , then If for every fixed x ∈ X , the mapping then is called time continuous. Time continuous semiflows are called 1-Lipschitz in finite horizon time if for every T > 0 , there exists a constant L = L T ≥ 0 such that for all x ∈ X , all t ∈ [0, T] , and all h ≥ 0 we have Next we will need bounds on the expansiveness of t . For this we introduce a modulus function.

(1) (0, x) = x for all x ∈ X, (2) (t + s, x) = (t, (s, x)) = (s, (t, x)) for all s, t ≥ 0 and all x ∈ X.

Definition 3 Let ∶ ℝ + ↦ ℝ + be a nondecreasing (growth) function satisfying (0) = 1 . We say that a semiflow is −equicontinuous on X if for every pair x, y ∈ X and all t ∈ ℝ + we have In particular, -equicontinuous semiflows are uniformly continuous in a finite horizon time. The function is said to be sub-multiplicative if for all u 1 , u 2 , … , u k ≥ 0 we have Clearly, if (t) = e t for some > 0 , then it satisfies (★ M ). It can be easily proved that condition (★ M ), with the extra assumption that is differentiable, reduces to the set of functions of the form Ce t .

Notice that -equicontinuity and the Lipschitz condition are related in the sense that as long as

The classical example of a semiflow is an ordinary differential equation

Example 2 Let us consider a simple ordinary differential equation

is an arbitrary positive, nonincreasing and continuous function. Its solutions (depending on the initial conditions x ∈ ℝ + ) constitute a semiflow t (x) = f (t, x) on ℝ + . Clearly the solutions satisfy

We can notice that

Our semiflow on the phase space satisfies condition ( ★ M ).

Example 3 The SIR model and its variations (SIS, SIRS, SEIR, SEIS Capasso 2008 , but see also Sects. 4 and 5 for details) is described by a system of ODEs. For example the SIS model has a general form S � (t) = f (t, S, I) , I � (t) = g(t, S, I) , assuming that for all t, S(t) + I(t) = 1 ≡ const. Thus, the SIS model can be considered as a semiflow t corresponding to those ODEs, where X = [0, 1] × [0, 1] ,

Similarly, the SIR model ODE system defines a semiflow on the space

• Take X = ℝ with the standard distance function | ⋅ | and define (t, x, v) = v t (x) = x + vt . Clearly, they form the classical affine flow. • Take X = ℝ 2 with the standard Euclidean distance and (t, (x, v), a) = a (t, (x, v)) = (a t 2 2 + vt + x, at + v) . This is another simple flow example.

We can see in all of the presented above examples that the transformations t depend on some parameters. Our extensions concern these parameters. Given a nonempty set of parameters, , we shall consider families of semiflows ℜ( ) = { ∶ ∈ } . In general, the space of parameters is endowed with a fixed -algebra G . It is assumed that all singletons are in the -algebra, { } ∈ G . Let us call ℝ + = [0, ∞) the time set, (X, ) the phase space, and ( , G) the control set. In this paper we shall operate mainly with finite subsets of a control set i.e.

Definition 4 Take, as mentioned above, * = {0, 1, … , M} to be a finite subset of .

with the additional restriction that n ≠ n+1 for all n. The family of all piecewise constant functions (the sequences of break points t 0 < t 1 < t 2 < ⋯ < t n → ∞ may vary with ) will be denoted by ℌ( * ) . For a fixed function ∈ ℌ( * ) , the family of transformations (where t ≥ 0)

is called a deterministically controlled (by ) semiflow. This could be considered to be a deterministic version of a renewal process. The class of all deterministically controlled by semiflows (i.e. transformations t (⋅), t ≥ 0 ) is denoted by ℜ ( * ) . We also define ℜ(ℌ( * )) = ⋃ ∈ℌ( * ) ℜ ( * ).

We now turn to discussing notations and other theoretical aspects. The proofs of the presented here lemmata and theorems can be found in the "Appendix".

the Dirac delta operation) becomes a metric space.

We skip the proof as it is obvious. Notice that the metric d r is equivalent to the L 1 -norm distance ∫ r 0 | (t) − (t)|dt . As the ranges of the functions from ℌ r ( * ) are confined to {0, 1, … , M} thus (on L 1 ([0, r]))

In particular, (ℌ r ( * ), d r ) is a separable metric space (not complete in general). We also notice that d r ( , ) = ({t ∈ [0, r] ∶ (t) ≠ (t)}).

Lemma 2 Let { t } t≥0, ∈ be the family of semiflows on (X, ) such that for finite * ⊆ the conditions (★) and ( ) hold, with a constant L and the growth function which is independent of ∈ * . Given , ∈ ℌ( * ) such that (0) = (0) , let w −1 = 0 < w 0 < w 1 < ⋯ be the sequence of real numbers, where the functions , stop being equal or alternatively start being equal

. For a fixed time parameter r > 0 let n * ∈ ℕ be such that w n * ≤ r < w n * +1 . Then for all x, y ∈ X we can bound if n * = 2j (the trivial case r ≤ w 0 will be discussed separately) and if n * = 2j − 1.

Before the reader turns to the proof in the "Appendix" we notice that if w 0 ≥ r , then evaluating ( r (x), r (y)) and using the bounds from Lemma 2 we will simply obtain ( r (x), r (y)) ≤ (r) (x, y) ; we remember that (0) = 1.

The next lemma is concerned with the case when (0) ≠ (0) and its proof is along the same lines as Lemma 2's. For the completeness of the paper we include almost all the details.

, ∈ be the family of semiflows on (X, ) such that, for some finite * ⊆ , the conditions (★) and ( ) hold, with a constant L and the growth function which is independent of ∈ * . Given , ∈ ℌ( * ) such that (0) ≠ (0) , let 0 = w 0 < w 1 < ⋯ be the sequence of real numbers, where the functions , stop being equal or alternatively start being equal (i.e.

. For a fixed time parameter r > 0 let n * ∈ ℕ be such that w n * ≤ r < w n * +1 . Then for all x, y ∈ X we have the bounds if n * = 2j and if n * = 2j − 1.

We mention that ( r (x), r (y)) ≤ 2Lr + (x, y) if w 1 ≥ r.

Replacing the condition (★) by the stronger one (★ M ), we obtain directly from the above lemmata:

, ∈ be the family of semiflows on (X, ) such that, for some finite * ⊆ , the conditions (★ M ) and ( ) hold, with a constant L and the growth function which is independent of ∈ * . Given , ∈ ℌ( * ) and a fixed time parameter r > 0 we have the bounds

In particular, the mapping ℌ r ( * ) ∋ ↦ r (x) ∈ X is (d r , ) continuous.

By ℌ r;m ( * ) we denote the set of all ∈ ℌ r ( * ) such that changes its values at most m times on the interval [0, r]. We obtain the next corollary.

Corollary 2 Let { t } t≥0, ∈ be the family of semiflows on (X, ) such that, for some finite * ⊆ , the conditions (★) and ( ) hold, with a constant L and the growth function which is independent of ∈ * . Given , ∈ ℌ r;m ( * ) and arbitrary x, y ∈ X we have the bounds and

The following theorem will be used in what will follow. It is a straightforward consequence of the already proved lemmata and corollaries.

Theorem 1 Let { t } t≥0, ∈ be the family of semiflows on (X, ) such that, for some finite * ⊆ , the conditions (★) and ( ) hold, with a constant L and the growth function which is independent of ∈ * . Then, for every fixed x ∈ X the mapping is (d r , ) continuous. Now let us relax and allow it to be the trajectories of a fixed continuous time homogeneous Markov chain (CTHMC) on * . It is well known (cf. Allen 2011, p. 200; Norris 2007 p. 73, 94) , that in this case we have • ( ) ∈ ℌ( * ) for almost all ∈ . Of course, the sequences of switching times t 0 ( ) < t 1 ( ) < t 2 ( ) < … are not fixed anymore (they do depend on a particular ∈ ), and therefore our control functions (⋅) = • ( ) may vary frequently. We describe this notion in all details by introducing first:

, which is a continuous time homogeneous Markov process on a measurable space ( * , G) . We denote Our definitions and notations concerning Markov processes are standard and are borrowed from Allen (2011), Gikhman and Skorokhod (1974) , Grimmett and Stirzaker (2009), Norris (2007) . Abbreviating, we shall denote the process { t } t∈[0,∞) by . In the case, when the control phase space * is finite, continuous time homogeneous Markov chains (CTHMC) are very simple and intuitive stochastic objects. Their dynamics and random evolution is fully governed by the so-called Q matrices Norris (2007) . A Q matrix is a real square matrix

such that all q j,k ≥ 0 , and for every j ∈ * we have ∑ M l=0,l≠j q j,l = q j,j . It follows that in the matrix Q the sum of elements in every row is 0. We will write q j instead of q j,j . The matrix Q is called the intensity matrix of the process as long as for every j ∈ * we have Formally, we may define the matrix Q as

We skip other formalisms and direct the reader to classical monographs (cf. Gikhman and Skorokhod 1974) , where all the essential results with clear proofs are presented. The intensity matrices (generators) are at the core of the analytical theory of Markov semigroups.

The trajectory of a process may be simulated using its Q matrix, by the classical Gillespie algorithm (Gillespie 1977) . Namely, a trajectory starting at t • = 0 , from an initial point 0 ∈ * , will stay in 0 for a random time 0 . Then, at the moment t 1 = 0 , it jumps with probability q 0 ,k q 0 to some k ≠ 0 . The scheme is repeated infinitely many times. As the process is Markovian, the distribution of , where ∈ * is exponential; i.e. P( > t) = e −q t , for t ≥ 0 (see Norris 2007, p. 94) . From the point of view of simulations, the full description of the trajectory (for a fixed ∈ ) is a sequence of pairs where n describe consecutive states of the trajectory of . The sequence of times t 0 < t 1 < t 2 < ⋯ represents instances of switching.

, t 2 = t 1 + 2 and so forth. We must notice that our presentation is a bit informal, as we skip over the technicality that all random variables (Markov moments) 0 , 1 , 2 , … are actually versions from the relevant families of independent copies. Obviously, t n ( ) is the time instance of the jump n ⇝ n+1 .

Given a CTHMC we shall denote by { n } n≥0 the embedded Markov chain and the sequence of the so-called inter-times of by { n } n≥0 . Of course n is a discrete time Markov chain on the control phase space * and n form a sequence of independent and nonnegative random variables, which describe the (waiting) time spent in n ( ) before jumping to n+1 ( ) . To keep formulae short (this repetition is consistent with the one introduced a while ago), let us denote t n ( ) = 0 ( ) + ⋯ + n ( ) , it is the time when the process t ( ) (after occupying

parameter n ( ) for the period n ( ) ) arrives at the state n+1 ( ) (being then "stuck" in it for the period n+1 ( )).

Let be a CTHMC defined on the probability space ( , F, P) with values in a finite control phase space ( * , G) and, as before, let ℜ( * ) = { (t, ⋅, ) ∶ t ≥ 0, ∈ * } be a family of semiflows on the phase space (X, ) . Then the mapping

is called a randomly controlled semiflow (which is a merged system of the family of the semiflows ℜ( * ) and the control CTHMC process = { t } t∈[0,∞) on the control phase space * ). We denote such a randomly controlled family of semiflows by ℜ( * ) ⋇ .

Each ∈ ℜ( * ) ⋇ is in fact a random process. For t ∈ ℝ + we introduce Of course ⋄+t • ( ) ∈ ℌ( * ) for P-almost all ∈ . It follows from our construction that (0, x, ⋅) = x . Moreover, if we denote t (x, ) = (t, x, ) then Theorem 2 Let { t } t≥0, ∈ be the family of semiflows on (X, ) such that, for finite * ⊆ , the conditions (★) and ( ) hold, with a constant L and the growth function which is independent of

} is a continuous time stochastic process on X with continuous trajectories.

We consider two simulation examples to illustrate the behaviour of the linear random dynamical system a t (x) = x + at . The first one is very simple and plays here an illustrative role. The second one is related to recently studied SARS-CoV-2 dynamics. One of the purposes of the simulation is to illustrate, similarly as will later be observed for SIR models, the limiting, as the frequency of change in the

and then inductively = (t − t n ( ), (t n ( ), x, ), n ( )) for t n ( ) ≤ t ≤ t n+1 ( ), ⋄+t r ( ) = r+t ( ) for all r ≥ 0.

driving process increases, behaviour of the considered a semiflow. Let a be a Markov process on a discrete state space {0, … , M} with transition rate matrix Q. Assume that a starts from 0. Let be the stationary measure of the Markov chain a. Define a * = lim

is a Markov process with transition rate matrix kQ, i.e. the matrix Q multiplied by the integer k. We assume that the parameter a is the realization of a Markov process on the state space {0, 1, 2} with rate matrix

We take this Q matrix for illustrative purposes, not due to some particular application. We first simulate the trajectory of a (k) for a large value of n = 100,000 . Then conditional on it, we calculate a t . Simulation of the Markov process is done by the Gillespie algorithm. We find a * by simulations. First from Q 1 we find the stationary distribution * as the first row of the matrix e 10,000 −1 , where is the diagonal matrix of eigenvalues of Q and is the matrix of eigenvectors (Thm. 21, Grimmett and Stirzaker 2009). The 10, 000 is the to represent ∞ , i.e. the limit of e tQ 1 as t → ∞ . From * we draw a large (10, 000) sample and take its average as a * . We simulate both a (k) and a * 100 times. We show the resulting collection of trajectories in Fig. 1 .

Next we aim at providing an example linked to the dynamics of the SARS-CoV-2 virus. Let the process a t represent a short term approximation of the growth Fig. 1 Simulated trajectories of a (k) t and a * t for Q 1 . The curves E a (k) t , a * t and E a * t essentially coincide with each other of the number of infected, i.e. at a stage of linear growth. The infected number, as described in the Introduction, is driven by a random factor. In the common SIR model (as will be described in Sect. 4) the number of infected is given by the differential equation and conditional on S(t) ≡ S we have the (conditional) solution where we can take S = R 0 . Hence, to compare with our model we would like a(t) = ( S − 1) , where some part of the right-hand side is time-dependent. Alternatively, we would like a * = ( S − 1) . To connect this with the dynamics of the SARS-CoV-2 virus one would need estimates of , R 0 for it and some population size S(t). Obtaining and R 0 from official statistics does not seem to be an easy task due to the nature of the collected data (Bartoszek et al. 2020) . In particular the number of confirmed cases strongly depends on the number of tests performed and furthermore one can suspect that the count will also depend on the testing strategy and type of test performed. However, one alternative is to use data collected from the cruise ship, Diamond Princess, which can be treated as a unique, naturally-occurring epidemiological study that could be useful for the prediction of parameters related to the behaviour of COVID-19. There were 3711 people onboard, 1045 crew (median age 36 years) and 2666 passengers (median age 69 year) (Moriarty et al. 2020 ). Of these 712 tested positive for the virus, 145 crew and 567 passengers 1 . The first case was observed at the end of January 2020 and by March 2020 all passengers and crew had disembarked (Moriarty et al. 2020) . Fourteen passengers died, with the first death on 10th February 2020 and last on 14th April 2020 (see Footnote 1). We need to have estimates of , the death rate, and R 0 . For the death rate we take (approximating the 64 days as 2 months), There are various estimates of R 0 for the Diamond Princess, 14.8, 1.78 ) and 2.28 (Zhang et al. 2020 ). Here we take the value of 3.7 (the reported by Liu et al. 2020 mean R 0 value from Wuhan and it was taken in Rocklöv et al. (2020) for a sensitivity analysis). We hence obtain a * DP = ((14∕567)∕2) ⋅ 3.7 ≈ 0.046 . Of course, in order to be able to carry out our simulation we need to have Q, of which a * is a function of. Let Q ⊂ ℝ 3×3 be the space of matrices that are constrained to have non-negative off-diagonals, each row summing to 0 and at least one non-zero off-diagonal entry in each row. Each Q ∈ Q will then be the transition rate matrix for some Markov process and have the a * Q parameter associated with. Our aim is to find a Q that minimizes

One cannot expect to have a unique Q , as multiple Qs can result in the same a * and hence we want to consider a whole set of viable Q s. We therefore take a Monte Carlo approach by running R's (R Core Team 2020) optim() function from 100 random seeds sampled from M . A given Q ∈ M is randomly chosen as follows. First for each row the position of the obligatory positive-off-diagonal entry is chosen. Then, from the remaining M − 1 (in this application only 1 entry remains, as M = 2 ) entries each one is chosen with probability 0.5 (we first draw from the binomial distribution the number of positive entries, and then using R's sample() draw the appropriate indices). The values of the chosen entries are drawn from the exponential with rate 1 distribution. Afterwords, the optimization of Eq. (3.1) is done using the Nelder-Mead method. We show the resulting collection of trajectories in Fig. 2 .

We finally turn to the investigation of trajectories of SIS (this Section) and SIR (Sect. 5) models controlled by a CTHMC with intensity matrix Q (and an initial distribution p) on finite state space {0, 1, … , M} (the numbers 0, 1, … , M here are just labels, the actual states of the process will be the values of the transmission rates). The classical SIS model for a constant population size Allen (2011), Kermack and McKendrick (1927) is described by a system of ODEs with constant

Qs optimized in order to have a * agreeing with a * DP coefficients (we remind the reader that without loss of generality we have assumed N = 1 , so in the more general case should be scaled by N and S(t) + I(t) = N):

The parameter is called the transmission rate-the number of contacts per time that results in an infection of a susceptible individual; is the recovery rate, where 1∕ is the average length of the infectious period.

As we have normalized out variables ( N = 1 ) we may consider solely the number of infected (and infectious) individuals I(t) at time t (as the number of susceptible individuals is directly given as S(t) = 1 − I(t)):

The above Bernoulli differential equation can be solved directly; the solution I(t) is given by The SIS epidemic model has an endemic equilibrium given by and it is positive if the basic reproduction number satisfies R 0 > 1 . In that case, the solutions approach the endemic equilibrium. If R 0 ≤ 1 solutions approach the infected-free state. The first SIS model is controlled by the CTHMC through the parameter (which, after further modifications, we denote by 1 ) which equals the value of the CTHMC.

For our simulation let us assume that the transmission rate = (t) is the realization of a CTHMC on the state space {0, 0.044, 0.11, 0.22} with intensity matrix From Q we find the stationary distribution ( * •Q = 0 ): * = (0.2, 0.4, 0.2, 0.2) . We note that the particular values of Q are chosen for illustrative purposes only. Figure 3 shows the trajectories of the piecewise-deterministic Markov process I(t)( ) (top panel), as well as the mean of L = 1000 trajectories (bottom panel),

superimposed on the ordinary deterministic solution I(t) with parameter corresponding to ≡ E[ ] = * •(0, 0.044, 0.11, 0.22). We postulate the limiting behavior of trajectories as we increase the intensities of jumps in CTHMC, i.e. the intensity matrix of the process, kQ changes as k → ∞ . Only numerical simulations were considered. However, our framework, considering limits of kQ as k → ∞ is consistent with known theorems describing conditions under which a sequence of Markov processes will converge to the solution of a system of first order ODEs (Kurtz 1971) . Figure 4 shows the trajectories of piecewise-deterministic Markov processes I(t)( ) with intensity matrices 10Q and 100Q superimposed on ordinary deterministic solution I(t) with parameter corresponding to * •(0, 0.044, 0.11, 0.22).

We now turn to considering a different stochastic modification of the SIS model. The spread of infection in the classical (deterministic) SIS model is associated with an incidence rate that is bilinear with respect to the number of susceptible Trajectories of the piecewise-deterministic Markov processes I(t)( ) for SIS models with random coefficients with intensity matrices 10Q and 100Q, and a comparison of the variance functions for models with intensity matrix Q, 10Q and 100Q and infected individuals. We propose a modification as we incorporate an "external" transmission of disease obtaining a non-multiplicative (affine) incidence rate (Korobeinikov and Maini 2005) . Motivations for such an external transmission have been provided in the Introduction and were also discussed in Sect. 3.

As a consequence we obtain the following ODE system, again for a constantsized population ( N = 1 ) SIS model:

where for all t > 0 , S(t) + I(t) = 1 ≡ const.

From the above relation we expand the equation for the number of infected individuals I(t):

We obtain again a first-order ODE, quadratic in the unknown function. However, our slight modification results in a significantly different "full-form" Riccati equation with constant coefficients (assuming 1 , 0 > 0).

We now turn to solving our ODE. We first notice that a particular solution i(t) is:

The "external" SIS epidemic model has an endemic equilibrium given by where = 1 − 0 − . Assuming 1 , 0 > 0 the endemic equilibrium is always positive, i.e. there is no infected free state. Solving by quadrature, substituting It can be proved that conditions ( ) and (★ M ) hold for the SIS model as the phase space is compact and (t ≥ 0 and x ∈ [0, 1] is an initial condition). Condition ( ) holds from the mentioned boundedness of derivatives. On the other hand, the condition (★ M ) can be derived from the SIS master equations (Eq. 4.1), and our assumption that S + I ≡ const(= 1) . Namely, we have Hence,

Similarly to the previous stochastic modification we implement the "external" rate 0 as a CTHMC with intensity matrix Q on the finite state space {0, 1, … , M} . Assume that the external incidence rate 0 = 0 (t) is the realization of a CTHMC on the state space {0, 0.044, 0.11, 0.01} with intensity matrix Q (as previously) and P( (0) = 0) = 1 . The other parameters for SIS model are as follows: I(0) = 0.01, 1 = 0.0836, = 0.05, K = 1. Figure 6 shows the trajectories of the process I(t)( ) superimposed on the ordinary deterministic solution I(t) with parameter 0 corresponding to * •(0, 0.044, 0.11, 0.01) . Again we find that with multiplied intensity matrices kQ we can observe the convergence to the deterministic 0 solution as k → ∞ , i.e. for high-frequency changes the in the underlying driving process, the observed process looks smooth.

From the formulation of the classical SIS model, assuming that at any given time t the number of infected individuals I(t) equals 0, we obtain a stable trivial equilibrium. However, it is not the case in our new model, as the positive "external" incidence rate makes it unstable. Figure 7 shows a sample trajectory presenting the origination of epidemics from a population free from infected individuals.

Similarly to the SIS (Sect. 4) models, we introduce an external stochastic infectious mechanism to the classical SIR model. We are unable to find analytical closed form solutions. The master equation of the SIR model guarantees that the generated by it semiflows allow for stochastic extensions. In order to obtain that the corresponding, controlled by a CTHMC, semiflows are stochastic processes, we simply explain that the crucial moment of the proof of our Thm 2 is the step where we show that the mapping Fig. 6 Trajectories of the piecewise-deterministic Markov processes I(t)( ) for SIS models with random coefficients with "external" incidence rate models with intensity matrices Q, 10Q, 100Q, red line denotes the deterministic solution of the ODE with parameter 0 is (d r , ) continuous (hence measurable). Now the continuity will follow from the boundedness of derivatives sup{|S � (t)|, |I � (t)|, |R � (t)|} < ∞ . Let us add that the corresponding vector fields on [0, 1] × [0, 1] × [0, 1] are furthermore smooth. Therefore, if the trajectories • ( ) are close in the metric d r on ℌ r;m then the final states (i.e. the ends of paths) • ( ) r (x) are also close. Let us stress that this section is strongly supported by numerics, instead of being a purely mathematical solution of the modified SIR model. As long as there is no explicit closed form solution it may be difficult to find (or estimate) the expansiveness of a semiflow. We can show that if x > y then |S t (x) − S t (y)| ≤ |xe − 0 t − ye −( 0 + 1 )t | . The expansiveness of I t (p) or R t (r) may be even more complicated.

Similarly to the SIS (Sect. 4) models, we introduce an external stochastic infectious mechanism to the classical SIR model. We again remind the reader that we have a normalized population size, N = 1.

where for all t > 0 , S(t) + I(t) + R(t) = 1 ≡ const. If 0 = 0 , then we obtain the classical SIR model. In the classical case it can be shown that lim t→∞ I(t) = 0 , and that the limits of S(t) and R(t) (as t → ∞ ) are finite but depend on the initial conditions. Assuming 0 = 0 , if the effective reproduction number Fig. 7 Trajectories of the piecewise-deterministic Markov processes I(t)( ) for for SIS models with random coefficients with "external" incidence rate models with intensity matrices Q, 0 (0) = 0 , 1 = 0.01 , = 0.05 , I(0) = 0 , K = 1 satisfies ≤ 1 , then there is no epidemic, and I(t)'s decrease to zero is monotonic. However, if > 1 , then I(t) first increases before decreasing to zero, a situation which we call an epidemic. Following the procedure used in Harko et al. (2014) we will now turn to analyzing the SIR model with external infection. We will not obtain a closed-form analytical solution, but reduce it to an Abel equation of the first kind. It follows from Eq. 

(5.5) (5.12)

(5.14) Polyanin and Zaitsev (1995) ). If 0 u = 0 the equation reduces to a Riccati-type equation, or if − 1 Du = 0 it reduces to a Bernoulli-type equation.

Equation (5.30) has no closed-form analytical solution. Introducing 0 > 0 to the SIR model results in a substantial change in our simulation process. We will rely on the R's (R Core Team 2020) package deSolve to solve the SIR ODE system, Eq. (5.1), between the interarrival times of the CTHMC.

(5.22) u u �� − (u � ) 2 + − 1 Du u u � + 0 u 2 = 0.

.

(u � ( (u))) 3 .

(5.26) u � ( (u(t))) = 1 (u(t)) ,

3 (u(t)) .

After introducing the "external" transmission rate 0 > 0 , we obtain only the trivial asymptotic equilibrium (S ⋆ , I ⋆ , R ⋆ ) = (0, 0, 1) . This is as all susceptible individuals will become eventually infected either via contact with the infected and infectious individuals (assuming I(t) > 0 ) or through the "external" transmission. Allowing for 0 to be described as a value of CTHMC provide us with a framework in which short "impulses" may switch the behavior of the solution (I(t), t > 0) between consecutive growth stages before decreasing to zero along the solution (S(t), t > 0) . Figure 8 shows a sample trajectories of SIR solutions assuming that 0 = 0 (t) is a Markov process with intensity matrix 0.1Q (for long inter-arrival times of the process) on a finite state space {0.001, 0.002, 0, 0.05} . We can clearly observe consecutive "waves" of the epidemic.

In line with our previous observations we can see that with the increase of the intensity matrix of the CTHMC kQ, as k → ∞ , the trajectories converge to the deterministic solution with parameter 0 corresponding to 0 ≡ E 0 = * •(0.001, 0.002, 0, 0.05) (see Fig. 9 ).

We end with a few historical remarks. The topic we presented is not new. The long-term behaviour of dynamical systems subjected to random perturbations was studied from different points of view (theoretical or practical computer simulations). It is quite impossible to give representative references, but we simply mention a few of them (see Benaïm and Strickler 2019; Chen-Charpentier and Stanescu 2010; Gray et al. 2012; Rami et al. 2014; Roberts 2017) . The common features of these models is that the dynamics is randomly interrupted so that the internal (very often hidden) state of the systems is changed, according to a law that depends on the state just before this intervention. In our work, the duration Describing the asymptotic behaviour of trajectories of randomly controlled semiflows seems to be a difficult task. There is no hope for achieving explicit analytical formulae, in particular, not even for deterministic SIR models, see Sect. 5. On the other hand, these notions provide theoretical models, which are recently a subject of interest for many research groups. In our work we employ computer techniques and find, by simulations, their statistical properties. The main goal is to find the dependence of these statistics for growing frequencies of the control processes (e.g. potential social contacts). This is reflected by considering sequences kQ of intensity matrices, where k → ∞ . For larger k the intensity matrix kQ will exhibit higher frequencies of change in the steering process. For special classes (in particular for SIR semiflows) the randomness of • ( ) t (x) = t (x, ) finally ( k → ∞ ) will be extinguished, approaching a deterministic semiflow for some specific mean ⋄ ∈ (not necessarily in * ).

For simplicity we considered control Markov processes that are irreducible on finite * . It follows that it has a unique stationary distribution * on * (see Allen (2011 ), Norris (2007 , p. 118). By the ergodic theorem (see Norris (2007) , p. 126) the trajectories t ( ) with probability 1 visit each fixed state ∈ * with frequency * ( ) . If the intensity coefficient k → ∞ , then on each (even very short) time interval (r, s) ⊂ [0, ∞) the control sequence (k) t ( ) switches parameters ∈ * accordingly to the vector * and many times (more and more when k is large), producing in some sense an averaged parameter ⋄ = ∑ ∈ * * ( ) . We notice that for * = {0, 1, 2, … , M} the expected value ⋄ generally does not belong to * any more. Therefore, the proposed approach is not universal, sometimes the parameters are not even numerical. In the considered SIS/SIR models, the parameters 0 (= • ( )) ∈ [0, M] behave in a linear (affine) Fig. 9 Trajectories of the piecewise-deterministic Markov processes I(t)( ) for SIR models with random coefficients with the "external" incidence rate SIR model with intensity matrix 10Q, 0 (0) = 0 , ▪ Proof of Lemma 3 If w 1 ≥ r (we apply conditions (★) and ( )), then Case I. Let n * = 2j . Reasoning as before we obtain ( w 1 (x), w 1 (y)) ≤ 2Lw 1 + (x, y) . Starting from t = w 1 we may apply the previous lemma as (w 1 ) = (w 1 ) . We simply replace r by r − w 1 and have Case II. n * = 2j − 1 . We begin as above ( r (x), r (y)) ≤ 2L((w 2j−1 − w 2j−2 ) (r − w 2j−1 ) + (w 2j−3 − w 2j−4 ) (r − w 2j−1 ) (w 2j−2 − w 2j−3 ) + (w 2j−5 − w 2j−6 ) (r − w 2j−1 ) (w 2j−2 − w 2j−3 ) (w 2j−4 − w 2j−5 ) + ⋯ + (w 1 − w 0 ) (r − w 2j−1 ) (w 2j−2 − w 2j−3 ) ⋯ (w 2 − w 1 )) + (r − w 2j−1 ) (w 2j−2 − w 2j−3 ) ⋯ (w 2 − w 1 ) (w 0 − 0) (x, y) .

( r (x), r (y)) ≤ ( r (x), 0 (x)) + ( 0 (x), 0 (y)) + ( 0 (y), r (y)) ≤ 2Lr + (x, y) .

( r (x), r (y)) = ( (⋅+w 1 )

+ (w 2j − w 2j−1 ) (w 2j−2 − w 2j−3 ) ⋯ (w 4 − w 3 ) (w 2 − w 1 )(2Lw 1 + (x, y)) = 2L((r − w 2j ) + (w 2j−1 − w 2j−2 ) (w 2j − w 2j−1 ) + (w 2j−3 − w 2j−4 ) (w 2j − w 2j−1 ) (w 2j−2 − w 2j−3 ) + ⋯ + (w 3 − w 2 ) (w 2j − w 2j−1 ) (w 2j−2 − w 2j−3 ) ⋯ (w 4 − w 3 ) + (w 1 − 0)) (w 2j − w 2j−1 ) (w 2j−2 − w 2j−3 ) ⋯ (w 4 − w 3 ) (w 2 − w 1 ) + (w 2j − w 2j−1 ) (w 2j−2 − w 2j−3 ) ⋯ (w 4 − w 3 ) (w 2 − w 1 ) (x, y) . 

We must show that for each fixed r ∈ [0, ∞) the mapping ∋ ↦ r (x, ) ∈ X is measurable. We illustrate this situation in Fig. 10 .

It follows from the construction that is a composition r (x, ) = • ( ) r (x) . If ∈ r;m , then the trajectory • ( ) , restricted to the interval [0, r], belongs to ℌ r;m ( * ) . By our Thm. 1, the mapping is (d r , ) continuous (so measurable). It remains to show that is measurable. Let A = {a 0 = 0, a 1 , a 2 , … } ⊂ [0, r] be a countable dense subset of [0, r] . Denote A k = {a 0 , a 1 , … , a k } . Consider a closed (for the relative d r topology) ball By separability, we can select a countable and dense set K = { 1 , 2 , … } in K r;m ( 0 , ) . We notice that Thus the trajectory • ( ) , restricted to the time domain [0, r], belongs to K r;m ( 0 , ) if and only if ( r (x), r (y)) = ( (⋅+w 1 ) r−w 1 ( w 1 (x)), (⋅+w 1 ) r−w 1 ( w 1 (y))) ≤ 2L((w 2j−1 − w 2j−2 ) (r − w 2j−1 ) + (w 2j−3 − w 2j−4 ) (r − w 2j−1 ) (w 2j−2 − w 2j−3 ) + ⋯ + (w 3 − w 2 ) (r − w 2j−1 ) (w 2j−2 − w 2j−3 ) ⋯ (w 4 − w 3 )) + (r − w 2j−1 ) (w 2j−2 − w 2j−3 ) ⋯ (w 4 − w 3 ) (w 2 − w 1 )(2Lw 1 + (x, y)) = 2L((w 2j−1 − w 2j−2 ) (r − w 2j−1 ) Clearly, { ∈ r;m ∶ a i ( ) = j (a i )} ∈ F ∩ r;m , as a i ∶ ↦ * is measurable. It follows that (see Fig. 10 ) r (x, ⋅) is F measurable. To end the proof we remark that the continuity of trajectories [0, ∞) ∋ t ↦ t (x, ) ∈ X is a trivial consequence of our construction of . ▪

An introduction to stochastic processes with applications to biology

A fractional model for the COVID-19 pandemic: Application to Italian data

Presumed asymptomatic carrier transmission of COVID-19

Are official confirmed cases and fatalities counts good enough to study the COVID-$19$ pandemic dynamics? A critical assessment through the case of Italy

Random switching between vector fields having a common zero

Mathematical structures of epidemic systems

Epidemic models with random coefficients. Mathematical and Computer Modelling

The theory of stochastic processes I

Exact stochastic simulation of coupled chemical reactions

The SIS epidemic model with Markovian switching

Probability and random processes

Exact analytical solutions of the susceptibleinfected-recovered (SIR) epidemic model and the SIR model with equal death and birth rates

Contribution to the mathematical theory of epidemics

Non-linear incidence and stability of infectious disease models

Limit theorems for sequences of jump Markov processes approximating ordinary differential processes

The reproductive number of COVID-19 is higher compared to SARS coronavirus

Markov chains

Handbook of exact solutions for ordinary differential equations

A language and environment for statistical computing. Vienna: R Foundation for Statistical Computing

Stability criteria for SIS epidemiological models under switching policies. Discrete and Continuous Dynamical Systems B

An epidemic model with noisy parameters

COVID-19 outbreak on the Diamond Princess cruise ship: Estimating the epidemic potential and effectiveness of public health countermeasures

Dynamical systems and population persistence

Stochastic dynamics II: Finite random dynamical systems, linear representation, and entropy production

Stochastic dynamics: Markov chains and random transformations. Discrete and Continuous Dynamical Systems

Estimation of the reproductive number of novel coronavirus (COVID-19) and the probable outbreak size on the Diamond Princess cruise ship: A data-driven analysis

We would like to thank two anonymous Reviewers whose comments significantly

fashion. Hence, one may evaluate the average ⋄ and ask whether the asymptotic behaviour (i.e. when k → ∞ ) of the stochastic dynamics is carried over.We used computer simulations to look at the high-frequency SIR properties. Namely, it looks like with growing frequencies the random trajectories become smoother and the total variance of the process decays to 0 and we seem to arrive at almost smooth graphs, characteristic for deterministic models. Hence, we can hypothesize that the high frequencies, i.e. short transition times, of the control process result in almost sure convergence of the controlled dynamical system.Let us finally comment that the high-frequency is imitated here by a direct increase of the intensity matrix, kQ. Of course, it would be desirable to work with more general approaches, e.g. considering kQ k , where for every k we only impose * •Q k = 0 . This is a direction for further theoretical work but here we dedicated to it a short, connected to the current spread of SARS-CoV-2 virus, Monte Carlo study. Interestingly, very recently a SIRD (SIR and Death) with random coefficients model was considered for describing the COVID-19 dynamics in Italy (Alòs et al. 2020 ). There they considered a fractional Brownian motion process for the infection parameter, in Eqs. (4.1) and (5.1). It is worth pointing out a key difference between modelling via a fractional Brownian motion and the random dynamical system approach presented here. A fractional Brownian motion in general will not have piecewise smooth (differentiable) trajectories. On the other hand, a random dynamical system, as described in this work, will.

Proof of Lemma 2 The functions , ∈ ℌ( * ) have finitely many jumps on the considered interval [0, r]. We recall that the total length of the intervals (w l , w l+1 ) on which the functions , differ is d r ( , ) (we are confined to the interval [0, r] only, thus the point w l+1 at the end of the final interval should be understood as r).Case I ( n * = 2j ). It follows from the condition (★) that as (t) = (t) for all t ∈ [0, w 0 ] (we remind again that if w 0 ≥ r , then ( r (x), r (y)) ≤ (r) (x, y)) ). On the interval [w 0 , w 1 ) the functions and differ, so applying condition ( ) we obtain (the constant L corresponds to T = r) Substituting 0 ← w 1 , w 0 ← w 2 , w 1 (x) ← x , and w 1 (y) ← y we obtain

Iterating this procedure we finally obtain This bound must be obviously confined to [0, r] , so we in fact can obtain Case II, n * = 2j − 1Taking t = r we obtain