On the exact reproduction number in SIS epidemic models with vertical transmission

Computational and Applied Mathematics (2023) 42:291 https://doi.org/10.1007/s40314-023-02424-5 On the exact reproduction number in SIS epidemic models with vertical transmission A. Gómez-Corral1·F. Palacios-Rodríguez2·M. T. Rodríguez-Bernal1 Received: 3 October 2022 / Revised: 12 June 2023 / Accepted: 1 August 2023 / Published online: 22 August 2023 © The Author(s) 2023 Abstract This paper proposes a bi-variate competition process to describe the spread of epidemics of SIS type through both horizontal and vertical transmission. The interest is in the exact reproduction number, Rexact,0, which is seen to be the stochastic version of the well-known basic reproduction number. We characterize the probability distribution function of Rexact,0 by decomposing this number into two random contributions allowing us to distinguish between infectious person-to-person contacts and infections of newborns with infective parents. Numerical examples are presented to illustrate our analytical results. Keywords Bi-variate competition process ·Reproduction number ·SIS epidemic model · Vertical transmission Mathematics Subject Classification 60J28 (Applications of continuous-time Markov processes on discrete state spaces) ·92D30 (Epidemiology) This work was supported by the Ministry of Science and Innovation (Government of Spain), under Project PGC2018-097704-B-I00. BA. Gómez-Corral [email protected] http://blogs.mat.ucm.es/agomez-corral/ F. Palacios-Rodríguez [email protected] M. T. Rodríguez-Bernal [email protected] 1Department of Statistics and Operations Research, Faculty of Mathematical Sciences, Complutense University of Madrid, Plaza de Ciencias 3, 28040 Madrid, Spain 2Department of Statistics and Operations Research, Faculty of Mathematics, University of Seville, Calle Tarfia s/n, 41012 Seville, Spain 123 291 Page 2 of 20 A. Gómez-Corral et al. 1 Introduction Vertical transmission of pathogens (Busenberg and Cooke 1992), such as viruses and bacteria, refers to generational transmission from colonized parents to their offspring. Several vertically transmitted infections are included in the term TORCH complex (see, e.g., Jaan and Rajnik 2021), which is related to congenital infections of Toxoplasmosis, Other infections (chickenpox caused by varicella zoster virus, chlamydia, hepatitis B, human immunodeficiency virus, parvovirus, syphilis, and Zika fever, among others), Rubella, Cytomegalovirus, and Herpes simplex. For some viruses—including human immunodeficiency viruses and hepatitis B virus—infections can occur by contact with or transfer of blood and body fluids (semen, vaginal fluids). This means that, in such cases, the generational transmission of the disease is possible at the same time that the person-to-person contact leads to the dissemination of the pathogen. Motivated by the above-mentioned pathogens, a variety of epidemic models with vertical and horizontal transmission have been analyzed, mainly under the deterministic perspective, with special emphasis on global stability properties such as the asymptotic stability of the disease-free equilibrium and the endemic equilibrium. Without claiming an exhaustive list of deterministic epidemic models involving vertical and horizontal transmission, we can mention a susceptible-infectious (SI) model (Kang and Castillo- Chavez 2014) of host–parasite interactions incorporating Allee effects and horizontal and/or vertical transmission, under the assumption that infectious individuals can experience pathogen-induced reductions in reproductive ability; a susceptible-infectious-susceptible (SIS) model with quarantined individuals and Beddington–DeAngelis incidence (Chen and Zhao 2020); a susceptible-exposed-infectious-removed (SEIR) model (Li et al. 2001); a susceptible-exposed-infectious-removed-susceptible (SEIRS) model (Gou and Wang 2007); and susceptible-infectious-removed (SIR) models with pulse vaccination (D’Onofrio 2005; Gao et al. 2007;Luetal.2002), and with nonlinear incidence rate, treatment, and vaccination for newborns of susceptible and recovered individuals (Hu et al. 2012), among others. The model of Naji and Hussien (2016) considers two strains with infections of SIRS type and of SIS type that spread through both horizontal and vertical transmission in the host population; see also Ref. Bichara et al. (2014), where the authors consider SIS and SIR models with n different pathogen strains and vertical transmission, including the possibility of infants with passive immunity in the framework of the SIR model. In the setting of plant–virus epidemiology, we refer the reader to the paper by Jeger et al. (1998), where an SEIR-type epidemic for the plant host is linked to an susceptible-exposed-infectious (SEI) model with vertical transmission for the insect vector population. The model of Jeger et al. (1998) describes the interaction of the virus with the vector and possible control options, including the removal and destruction of diseased plants and reduction of the vector-population size, for example, by insecticide or vegetation management. A non-linear model is introduced by Naresh et al. (2006) to study the transmission of HIV/AIDS into a population with varying size with vertical transmission and other demographic and epidemiological aspects. In an attempt to investigate how the unpredictability of person-to-person contact affects the spread of a disease, Kiouach and Sabbar (2018) assume that the contact rate is perturbed by Gaussian white noise and establish sufficient conditions for persistence and extinction of the disease in a stochastic SIRS model with vertical transmission and transfer from infectious to susceptible compartments. In the setting of the SIR model with vertical transmission and media coverage, stochastic perturbations are introduced by Wang et al. (2020), who prove existence and uniqueness of the positive solution, and extinction and persistence in mean. 123 On the exact reproduction number in SIS epidemic models... Page 3 of 20 291 In analyzing random perturbations in the SIS model with vertical transmission and immigration of susceptible individuals, and emigration of susceptible and infectious individuals, Zhang et al. (2018) use suitably defined Lyapunov functions to derive an ergodic stationary distribution, and conditions for the extinction and persistence of the pathogen. In this paper, we focus on a Markov chain model to describe the spread of epidemics of SIS type through both horizontal and vertical transmission. Accordingly, the underlying Markov chain is seen as a competition process (Iglehart 1964;Reuter1961) on the positive orthant N0×N0, where the point (i,s)represents the size of the subpopulations of infectious and susceptible individuals, and the one-step transitions from (i,s)are allowed only to a neighboring point (i,s)∈{(i−1,s), (i−1,s+1), (i,s−1), (i,s+1), (i+1,s− 1), (i+1,s)}; as a result, this process can be thought of as the natural generalization of the uni-variate birth–death process arising from the standard SIS epidemic model. Specifically, the interest is in characterizing the probabilistic behavior of the exact reproduction number Rexact,0(Artalejo and Lopez-Herrero 2013; Economou et al. 2015, Section 3.3), which is decomposed here into two random contributions that record infectious contacts between individuals and infections of newborns with infectious parents. Our methodology is based on the use of finite quasi-birth–death (QBD) processes and related algorithms to evaluate their stationary vector (Akar et al. 2000; De Nitto Personè and Grassi 1996; Gaver et al. 1984, Section 2), first-passage times and sojourn times (Gaver et al. 1984, Sections 3–4; Gómez-Corral et al. 2020), and perturbation properties (Gómez-Corral and López-García 2018). In the deterministic setting, the basic reproduction number R0is by far the most widely used index in mathematical epidemiology (Li et al. 2011; Roberts 2007), mainly due to the simplicity of its definition and the fact that it is usually a threshold value establishing that a disease persists only if R0>1 (van den Driessche and Watmough 2002). It is also important to remark its role in the validity of the final size equation and the uniqueness of its solution (Ma and Earn 2006), and the control effort required to eradicate the pathogen, among other aspects. When a stochastic approach is adopted, it is well known (Artalejo and Lopez-Herrero 2013) that the value of R0overestimates the reproductive potential of a pathogen, especially when the involved community has a small size. An appropriate alternative to R0—in particular, for Markov chain models in epidemics—is the exact reproduction number Rexact,0,which is introduced in Artalejo and Lopez-Herrero (2013) by translating the definition of R0into a random variable. Among the most relevant features of Rexact,0, compared to R0,arethat Rexact,0eliminates the effect of repeated infectious contacts; it is not necessarily defined at the time of invasion, but at any later time; and it provides a more accurate index to quantify the transmission of a disease. However, the analytical treatment of Rexact,0requires knowing apriorithe underlying parameters of the epidemic model. In computing the mass function of Rexact,0and its moments, the most intensive effort lies in the inversion of matrices and related memory requirements, whence the numerical solution of Rexact,0may experience instability with a large population size. These features and limitations show the adequacy of Rexact,0to quantify, in a more realistic way than R0, the transmission potential of a disease in small communities sharing confined spaces, such as families and intensive care units (Economou et al. 2015), a mini-community design to assert the indirect effect of vaccination (Fernández-Peralta and Gómez-Corral 2021), and hospital ward (Chalub et al. 2023), among others. The organization of the paper is as follows. In Sect. 2, we first give the mathematical description of the SIS epidemic model with vertical and horizontal transmission, and then introduce a bi-variate competition process to describe the number of susceptible and infectious individuals at time t. In Sect. 3, we characterize the joint probability law of the exact 123 291 Page 4 of 20 A. Gómez-Corral et al. Fig. 1 Diagram of transitions between compartments in the SIS epidemic model with vertical and horizontal transmission number RH exact,0of secondary infections by contact between susceptible individuals and a (predetermined) initially infectious individual, and the number RV exact,0of infectious newborns with the initially infectious individual as parent. Based on the tridiagonal-by-blocks form of the underlying q-matrix, we derive an iterative procedure for computing the joint probability law of (RV exact,0,RH exact,0)by using first-step analysis and block-Gaussian elimination. In Sect.4, our analytical results are exemplified with some numerical experiments. A brief discussion is finally presented in Sect.5. 2 The mathematical model description Consider a closed and homogeneously well-mixed population of N(t)individuals that can be partitioned into two compartments or subpopulations of susceptible (S), and infectious (I) individuals, with densities, respectively, denoted by S(t)and I(t)at time t;seeFig.1.A susceptible individual is assumed to be healthy, but can become infectious by contact with an infectious individual according to a Poisson process of rate μ>0 in such a way that it should remain infectious during an exponentially distributed time of parameter γN>0. The natural recovery period (with mean length γ−1 N) of an infectious individual can be effectively shortened due to health-care actions with treatment function η(i)=γTmin{i,C},where γT≥0 denotes the treatment rate,1and Cis the threshold at which the health-care system reaches its maximal capacity. This means that treatment increases linearly with the number i of infectious individuals before reaching the maximum resources and is constant afterward. The model incorporates the birth and the death of individuals, and the vertical transmission of the disease. Therefore, the total population size N(t)=S(t)+I(t)is not constant, but it is assumed that N(t)∈{0,1, ..., M}in the sense that limited resources in the ecosystem prevent population growth above a certain maximum size M<∞. An infectious individual gives birth to either a susceptible individual with rate pβIor an infectious one with rate (1−p)βI, where βI>0 is the birth rate and p∈[0,1]represents the proportion of the offspring of infectious parents that are susceptible individuals. Only susceptible individuals are assumed to be born with susceptible parents, which occurs at rate βS>0. There is a natural death rate αS>0 for susceptible individuals, and infectious individuals face death due to either natural causes or the disease with joint death rate αI>0. The underlying processes governing infectious contacts and infections of newborns, recovery times, and birth and death events are assumed to be mutually independent. 1The selection γT=0 amounts to a treatment that alleviates the severe symptoms, but does not reduce the infection time. 123 On the exact reproduction number in SIS epidemic models... Page 5 of 20 291 Fig. 2 Transitions from any state (i,s)∈S\l(0) The birth rates βSand βIare not necessarily assumed equal, which will allow us to analyze the influence of birth control measures, for example, on the subpopulation of infectious individuals, taking decreasing values of βI, for a fixed value of βS. The choice αI=α I+α I allows us to distinguish between the death of infectious individuals due to natural causes or due to the disease, with αS=α Iif the death of an individual does not depend on the compartment and, consequently, αI>α Swhen the disease can cause the death of infectious individuals. We note that the model coincides with the SIS epidemic model with demography if p=1, which results in a uni-variate birth–death process on the state space {0, ..., M}. The standard SIS epidemic model (Allen 2003, Section 7.3.2) is derived when the births and deaths are ignored with the selection αS=αI=βS=βI=0. The above Markovian formulation yields a bi-variate competition process X= {(I(t), S(t)) :t≥0}(Reuter 1961) taking values in the finite reducible state space S={(i,s)∈N0×N0:i+s∈{0, ..., M}}, whose non-vanishing infinitesimal transition rates (Fig.2) from state (i,s)to state (i,s)are given by q(i,s),(i,s)= ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ αIi,if (i,s)=(i−1,s), αSs,if (i,s)=(i,s−1), γNi+η(i), if (i,s)=(i−1,s+1), μis,if (i,s)=(i+1,s−1), (1−p)βIi,if (i,s)=(i+1,s), βSs+pβIi,if (i,s)=(i,s+1), for states (i,s), (i,s)∈Swith (i,s)= (i,s). These rates are related, respectively, to the death of an infectious individual, the death of a susceptible individual, the recovery of an infectious individual, an infectious contact between the compartments of susceptible and infectious individuals, the birth of an infectious individual, and the birth of a susceptible individual. We also consider the value q(i,s),(i,s)=−q(i,s),whereq(i,s)is the rate of the sojourn time of process Xin state (i,s);thatis q(i,s)=(αI+(1−δi+s,M)βI+γN)i+(αS+(1−δi+s,M)βS)s+η(i)+μis, 123 291 Page 6 of 20 A. Gómez-Corral et al. for (i,j)∈S,whereδa,brepresents the Kronecker delta. Note that the state (0,0)is absorbing, and Xevolves as a birth–death process on the class {(0,s):s∈{0, ..., M}} of infection-free states, with birth parameters {βSs:s∈{1, ..., M−1}} and death parameters {αSs:s∈{1, ..., M}}. For later use, we use a lexicographical ordering of states and decompose the state space Sby levels into ∪M i=0l(i),wheretheith level is defined as the subset l(i)={(i,s):s∈{0, ..., M−i}}, for i∈{0, ..., M}. This labeling of states allows to reformulate the process Xas a finite QBD process (Gaver et al. 1984;Gómez-Corraletal.2020) with q-matrix Q= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ Q0,0Q0,1 Q1,0Q1,1Q1,2 ......... QM−1,M−2QM−1,M−1QM−1,M QM,M−1QM,M ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ,(1) where the sub-matrices Qi,iare seen to be of dimension (M−i+1)×(M−i+1), for integers i∈{max{0,i−1},i,min{i+1,M}} and i∈{0, ..., M}. In particular, the non-vanishing entries of Qi,i−1are related to the removal of an infectious individual due to its recovery or death; the diagonal entries of Qi,iare given by −q(i,s)for states (i,s) with s∈{0, ..., M−i}, whereas the non-vanishing off-diagonal entries of Qi,iare linked to the birth and the death of a susceptible individual; and the non-vanishing entries of Qi,i+1 are associated with new infections due to the horizontal and the vertical transmission of the disease. Expressions for these sub-matrices are summarized in Appendix A. An immediate consequence of the resource-limited condition (i.e., N(t)∈{0,1, ..., M} with M<∞) is the fact that the length of an outbreak is a non-defective random variable. This means that, from any initial state (i,s)∈S, process Xwill reach the class of infection-free states almost surely. Indeed, the mean time to the first entrance into the class of infection-free states is finite, regardless to the initial state (i,s)∈S. Remark 1 The dynamics of process Xin the case M=∞is related to a QBD process with infinitely many states, for which further work must be done to study the occurrence of the first passage to the class {(0,s):s∈N0}of infection-free states. While it is outside the scope of this work, we briefly point out here that, by Gómez-Corral et al. (2023, Theorem 1), the resulting QBD process with infinitely many states is always regular, whence inequalities max{βI,βS}<min{αI,αS},andmin{βI,βS}≤αI+αSif max{βI,βS}=min{αI,αS}, can be seen as sufficient for the first passage to the subset {(0,s):s∈N0}of disease-free states to occur with certainty, and in a finite mean time, from any initial state (i,s)∈N0×N0; see, e.g., Gómez-Corral et al. (2023, Theorem 3) and Reuter (1961, Theorems 2-3). 3 An exact measure of the disease spread Let us consider an invasion time (i.e., I(0)=1andS(0)=s0with s0∈{0, ..., M−1})and define the exact number Rexact,0of secondary infections generated by the initially infectious individual before it recovers or dies (Artalejo and Lopez-Herrero 2013), which should occur at a certain random time τ. Specifically, time τis an exponentially distributed random variable with mean value (αI+γN+(1−δC,0)γT)−1. The most remarkable feature of process Xis 123 On the exact reproduction number in SIS epidemic models... Page 7 of 20 291 Table 1 Prevalence estimates of congenitally acquired disease, measured in terms of cases per 100 live births, for some TORCH diseases Disease 1 −pSource References Toxoplasmosis 0.01–0.033% United States Neu et al. (2015) HIV 31–36% 2005, World Neu et al. (2015) Treponema pall. (Syphilis) 0.0078% 2012, United States Neu et al. (2015) Rubella 0.0004–0.002% 1993–1997, Australia Sullivan et al. (1999) Cytomegalovirus 0.08% United States Neu et al. (2015) 0.6–0.7% World Marsico and Kimberlin (2017) Herpes simplex 0.0103% 2010–2015, World Bhatta et al. (2018) 0.0047% 2006–2011, Netherlands Bhatta et al. (2018) the vertical and the horizontal transmission of the disease, whence we decompose Rexact,0 into two random contributions RV exact,0and RH exact,0according to the fact that infections are linked, respectively, to infectious newborns with the initially infectious individual as parent and to infectious contacts between the compartment of susceptible individuals and the initially infectious individual, occurring both events before time τ. To prevent transmission of the disease at the initiation of an outbreak, the decomposition of Rexact,0into RV exact,0 and RH exact,0is suited to prioritize birth control measures, especially from infectious parents, over contingency measures against person-to-person contact—e.g., social distancing, and handwashing—or vice versa, and to analyze the subsequent consequences that limitations on these prophylactic interventions could cause on the relationships between individuals and even the survival of the population. For instance, when the magnitude of Rexact,0—in terms of sample paths—or its expectation is greater than one, the values of RV exact,0, and of RH exact,0 will make it possible to quantify, in a differentiated way, the effects of specific interventions on the dynamics of the population. The probability 1 −pof infectious offspring (Table 1) is closely related to the disease, whence no action can be taken on the value of 1 −pin the framework of the above-mentioned interventions. Let I∗denote the status of the initially infectious individual at time τ+0 (i.e., immediately after time τ). In particular, I∗=Dif the initially infectious individual dies before recovering; I∗=Nif it recovers due to natural causes; and I∗=Tif it recovers due to treatment. The aim is to characterize the joint probability law of (RV exact,0,RH exact,0)by evaluating the conditional probabilities P(1,s0)(i∗,rV,rH)=PI∗=i∗,RV exact,0=rV,RH exact,0=rH  (I(0), S(0)) =(1,s0), for i∗∈{D,N,T}and pairs (rV,rH)∈N0×N0. It is clear that RV exact,0and RH exact,0are inherently linked to events taking place during the time interval [0,τ) and, consequently, the conditional probabilities P(1,s0)(i∗,rV,rH),for status i∗∈{D,N,T}and pairs (rV,rH)∈N0×N0, satisfy the equalities αI αI+γN+(1−δC,0)γT = (rV,rH)∈N0×N0 P(1,s0)(D,rV,rH), (2) γN αI+γN+(1−δC,0)γT = (rV,rH)∈N0×N0 P(1,s0)(N,rV,rH), (3) 123 291 Page 8 of 20 A. Gómez-Corral et al. Fig. 3 An iterative procedure to compute the family of conditional probabilities {P(i,s)(i∗,rV,rH):(i,s)∈ S\l(0)}, for a predetermined status i∗and a fixed pair (rV,rH), in the special case rV=rH+2 (1−δC,0)γT αI+γN+(1−δC,0)γT = (rV,rH)∈N0×N0 P(1,s0)(T,rV,rH). (4) These equalities are derived by observing that the left-hand side of (2)–(4) corresponds to the marginal mass function {P(I∗=i∗|(I(0), S(0)) =(1,s0)) :i∗∈{D,N,T}}. In evaluating the probability P(1,s0)(i∗,rV,rH), for a predetermined status i∗∈{D,N,T} and a fixed pair (rV,rH)∈N0×N0,wefirstdefine,forstates(i,s)∈S\l(0)and any time t< τ, the conditional probability P(i,s)(i∗,r V,r H)that the initially infectious individual gives birth to r Vinfectious offspring and generates r Hinfectious contacts with susceptibles during the residual interval [t,τ), and its status at time τ+0isi∗, given that (i,s)is the current state of Xand the infectious period of the initially infectious individual is in process at time t.We then derive an iterative procedure that, starting from the families of conditional probabilities {P(i,s)(i∗,0,0):(i,s)∈S\l(0)},{P(i,s)(i∗,r V,0):r V∈{1, ..., rV},(i,s)∈S\l(0)}and {P(i,s)(i∗,0,r H):r H∈{1, ..., rH},(i,s)∈S\l(0)}, evaluates the conditional probabilities in {P(i,s)(i∗,rV,rH):(i,s)∈S\l(0)}in terms of {P(i,s)(i∗,rV−1,rH):(i,s)∈S\l(0)} and {P(i,s)(i∗,rV,rH−1):(i,s)∈S\l(0)}, for every status i∗∈{D,N,T}and pair (rV,rH)∈N×N;seeFig.3. 3.1 A first-step analysis of process X We can derive a system of linear equations for the conditional probabilities P(i,s)(i∗,rV,rH), for (i,s)∈S\l(0),i∗∈{D,N,T}and (rV,rH)∈N0×N0, using the first-step analysis. To be concrete, it is found that P(i,s)(D,0,0)=αI q(i,s) +(1−δi,1)αI(i−1) q(i,s) P(i−1,s)(D,0,0) +(1−δi,1)γN(i−1)+(1−δC,0)γTmin{i−1,C−1} q(i,s) P(i−1,s+1)(D,0,0) +(1−δs,0)αSs q(i,s) P(i,s−1)(D,0,0) 123 On the exact reproduction number in SIS epidemic models... Page 9 of 20 291 +(1−δi,M)(1−δs,M−i)βSs+pβIi q(i,s) P(i,s+1)(D,0,0) +(1−δi,1)(1−δs,0)μ(i−1)s q(i,s) P(i+1,s−1)(D,0,0) +(1−δi,1)(1−δi,M)(1−p)βI(i−1) q(i,s) P(i+1,s)(D,0,0), (5) for states (i,j)∈S\l(0). Similarly, the conditional probabilities P(i,s)(N,0,0)and P(i,s)(T,0,0),for(i,j)∈S\l(0), are seen to satisfy (5) with the term q−1 (i,s)αIreplaced by q−1 (i,s)γNand q−1 (i,s)(1−δC,0)γT, respectively. For any status i∗∈{D,N,T}and pairs (rV,rH)∈N0×N0\{(0,0)}, it is also seen that P(i,s)(i∗,rV,rH)=(1−δi,1)αI(i−1) q(i,s) P(i−1,s)(i∗,rV,rH) +(1−δi,1)γN(i−1)+(1−δC,0)γTmin{i−1,C−1} q(i,s) P(i−1,s+1)(i∗,rV,rH) +(1−δs,0)αSs q(i,s) P(i,s−1)(i∗,rV,rH) +(1−δi,M)(1−δs,M−i)βSs+pβIi q(i,s) P(i,s+1)(i∗,rV,rH) +(1−δi,M)(1−δs,0)μs q(i,s) P(i+1,s−1)(i∗,rV,rH−1) +(1−δi,1)(1−δs,0)μ(i−1)s q(i,s) P(i+1,s−1)(i∗,rV,rH) +(1−δi,M)(1−p)βI q(i,s) P(i+1,s)(i∗,rV−1,rH) +(1−δi,1)(1−δi,M)(1−p)βI(i−1) q(i,s) P(i+1,s)(i∗,rV,rH), (6) for states (i,j)∈S\l(0). Equations (5)and(6) can be easily written in matrix form using column vectors pi(i∗,rV,rH),fori∈{1, ..., M}, with (1+s)th entry P(i,s)(i∗,rV,rH)for integers s∈{0, ..., M−i}, and a suitable decomposition of sub-matrices Qi,i−1and Qi,i+1. Specifically, the rate sub-matrix Qi,i−1is written in terms of Qi,i−1=QD i,i−1+QN i,i−1+QT i,i−1+Q∗ i,i−1,i∈{1, ..., M}, where QD i,i−1,QN i,i−1,andQT i,i−1consist of the infinitesimal transition rates from states (i,s) with s∈{0, ..., M−i}associated, respectively, with the death, the natural recovery, and the recovery due to treatment of the initially infectious individual, and Q∗ i,i−1is related to transitions associated with the death and recovery of the remaining i−1 infectious individuals. The rate sub-matrix Qi,i+1is written in a similar manner as follows: Qi,i+1=QV i,i+1+QH i,i+1+Q∗ i,i+1,i∈{1, ..., M−1}, where QV i,i+1and QH i,i+1record the infinitesimal transition rates associated, respectively, with the vertical and the horizontal secondary infections generated by the initially infectious individual before time τ,andQ∗ i,i+1records the infinitesimal transition rates associated with secondary infections generated by the remaining i−1 infectious individuals. 123 291 Page 16 of 20 A. Gómez-Corral et al. Fig. 8 Diagram of transitions between compartments in an SIS epidemic model with vertical and horizontal transmission, and pre- and post-exposure vaccination sider here the possibility of enriching the model by allowing vaccination (Fig.8), and we exemplify how the QBD formulation is the key in the derivation of the joint distribution of (RV exact,0,RH exact,0). A third compartment of vaccinated (V) individuals must be incorporated into the model, and the infinitesimal dynamics of the resulting three-variate competition process X={(I(t), S(t), V(t)) :t≥0}(Iglehart 1964) are governed by the following non-vanishing transition rates from state (i,s,v)to state (i,s,v): q(i,s,v),(i,s,v)= ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ αIi,if (i,s,v)=(i−1,s,v), αSs,if (i,s,v)=(i,s−1,v), αVv, if (i,s,v)=(i,s,v−1), (1−¯p)(γNi+η(i)), if (i,s,v)=(i−1,s+1,v), ¯p(γNi+η(i)), if (i,s,v)=(i−1,s,v+1), λs,if (i,s,v)=(i,s−1,v+1), μis,if (i,s,v)=(i+1,s−1,v), (1−p)βIi,if (i,s,v)=(i+1,s,v), ¯p(βSs+βVv) +¯ppβIi,if (i,s,v)=(i,s,v+1), (1−¯p)(βSs+βVv) +(1−¯p)pβIi,if (i,s,v)=(i,s+1,v), for states (i,s,v),(i,s,v)∈Swith (i,s,v)= (i,s,v) and S={(i,s,v) ∈N0× N0×N0:i+s+v∈{0, ..., M}},andq(i,s,v),(i,s,v) =−q(i,s,v),where q(i,s,v) =(αI+(1−δi+s+v,M)βI+γN)i+(αS+(1−δi+s+v,M)βS+λ)s +(αV+(1−δi+s+v,M)βV)v +η(i)+μis, for (i,s,v) ∈S. In these expressions, αVis the death rate of vaccinated individuals, βV is the birth rate of healthy individuals from vaccinated parents, and λis the rate that a susceptible individual becomes vaccinated by administration of a therapeutic (post-exposure) vaccine. Probabilities ¯pand ¯p, with ¯p,¯p∈[0,1], are linked to prophylactic (pre-exposure) vaccination interventions at birth of healthy individuals from susceptible and vaccinated parents, and infectious parents, respectively, and probability ¯p, with ¯p ∈[0,1],is related to the use of a prophylactic vaccine after recovery of infectious individuals. Vaccines are assumed to be perfect; i.e., vaccines are fully effective and vaccinated individuals acquire lifelong immunity against the pathogen. As is the case of process X,theq-matrix of Xis seen to have the structured form of Qin Eq. (1) using the labeling of states 123 On the exact reproduction number in SIS epidemic models... Page 17 of 20 291 S=∪ M i=0l(i)with levels l(i)={(i,s,v) :s+v∈{0, ..., M−i}},fori∈{0, ..., M}. This means that, for process X, sub-matrices Qi,i−1and Qi,i+1in (1) record infinitesimal rates associated with the recovery and the death of infectious individuals, and new infections, respectively, and diagonal entries of Qi,iare specified by −q(i,s,v),for(i,s,v) with s+v∈{0, ..., M−i}. Non-vanishing off-diagonal entries of Qi,icorrespond to the birth and the death of susceptible and vaccinated individuals, and the therapeutic use of a vaccine on susceptible individuals. The column vectors pi(i∗,rV,rH)for process X, for integers i∈{1, ..., M}and i∗∈{D,N,T}, are defined as pi(i∗,rV,rH)in Section 3.1 with the probabilities P(i,s)(i∗,rV,rH)replaced by the conditional probabilities P(I∗=i∗,RV exact,0=rV,RH exact,0=rH|(I(0), S(0), V(0)) =(i,s,r)), for states (i,s,v)∈S\l(0). Equations (9)and(10) are then found to determine the joint probability law of (RV exact,0,RH exact,0)by adapting our approach in Sections 3.1 and 3.2 to process X;in particular, matrices Qi,i−1(respectively, Qi,i+1)forXmust be appropriately decomposed by distinguishing, in a similar manner to X, between the death, the natural recovery, and the recovery due to treatment of the initially infectious individual (respectively, the vertical and the horizontal secondary infections generated by the initially infectious individual). The use of any prophylactic or therapeutic vaccine will decrease the opportunities of the initially infectious individual to infect healthy individuals, whence model Xis expected to predict a lower transmission of the pathogen at an early state of the disease than model X; i.e., events {RV exact,0∈{0, ..., r V}} and {RH exact,0∈{0, ..., r H}}, for certain relatively small values r Vand r H, are expected to occur more likely in process Xthan in process Xand, on at least first observation, the tail of the distribution of RV exact,0(respectively, RH exact,0)for process Xis expected to be heavier than that of RV exact,0(respectively, RH exact,0) for process X. Further work must be done to predict the effects of a prophylactic intervention—defined by the selection λ=0 and at least one strictly positive value of ¯p,¯pand ¯p—compared to a therapeutic intervention—defined by a strictly positive value of λand ¯p=¯p=¯p =0— especially to derive certain thresholds for the underlying probabilities ¯p,¯pand ¯p,andthe rate λthat make the use of a prophylactic vaccine preferable over the use of a therapeutic one, and vice versa. In the Markov chain model, person-to-person contacts are governed by mutually independent Poisson processes of common rate μ, regardless of whether they result in a new infection or not. This means that, since the Poisson process has independent and stationary increments (Kulkarni 2020, Chapter 5), individuals are assumed to together constitute a homogeneously well-mixed population with varying size over time. Another aspect of interest is how to preserve the Markovian formulation of the model when the exponential regime—governing infectious contacts, recovery times, and birth and death events—is generalized to non-exponential distributional assumptions. For instance, a Markov-modulated version of process Xcan be defined by replacing the exponentially distributed recovery times with mean γ−1 Nby phase-type random variables (He 2014, Chapter 1) of order m. Instead of states (i,s)for X, this results in states (i,s;y)for an augmented process (X,Y)with vector y=(y1, ..., ym)—with m j=1yj=i—recording the number yjof phase-type recovery times at phase j,for j∈{1, ..., m}, at an arbitrary time; note that other Markov-modulated versions can be obtained from the replacement of exponentially distributed times until the occurrence of individuals’ birth and death by phase-type random variables, and/or the Poisson processes for infectious contacts by Markovian arrival streams (He 2014, Chapter 2). Although process (X,Y)can be thought of as a finite QBD process and the approach in Section 3 can be readily adapted, the memory requirements are so demanding that the numerical solution for the mass function of (RV exact,0,RH exact,0)will become rapidly prone to numerical 123 291 Page 18 of 20 A. Gómez-Corral et al. instability due to the cardinality of the underlying state space, as observed in Ref. Almaraz and Gómez-Corral (2018) for the simpler SIR model. In analyzing the exact reproduction number in SIS epidemic models with vertical transmission and a general infection period distribution, the theory on piecewise-deterministic Markov processes will be requested in a similar manner to the study of secondary cases in the framework of the SIR model (Gómez- Corral and López-García 2017). Acknowledgements The authors would like to thank the editor and four anonymous referees for their constructive comments, which have improved the presentation of the paper. Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. Data Availability No data was used for the research described in the article. Details on numerical experiments in Sect. 4will be made available on request. Declarations Conflict of interest The authors declare no potential conflicts of interest. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Appendix A: The q-matrix of process X We construct the structured q-matrix in Eq. (1) using a lexicographical labeling of states. This yields the following expressions for sub-matrices Qi,i,fori∈{max{0,i−1},i,min{i+ 1,M}}: (i) For i=0, the sub-matrix Q0,0has the form Q0,0= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 αS−q(1,0)βS 2αS−q(2,0)2βS ......... (M−1)αS−q(M−1,0)(M−1)βS MαS−q(M,0) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , and, since l(0)is a closed class of states, Q0,1=0(1+M)×M,where0 a×bdenotes the null matrix of dimension a×b. (ii) For i∈{1, ..., M}, the sub-matrices Qi,1−1and Qi,iare given by Qi,i−1= ⎛ ⎜ ⎜ ⎜ ⎝ αIiγNi+η(i) αIiγNi+η(i) ...... αIiγNi+η(i) ⎞ ⎟ ⎟ ⎟ ⎠ , 123 On the exact reproduction number in SIS epidemic models... Page 19 of 20 291 Qi,i= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ −q(i,0)pβIi αS−q(i,1)βS+pβIi 2αS−q(i,2)2βS+pβIi ......... (M−i−1)αS−q(i,M−i−1)(M−i−1)βS+pβIi (M−i)αS−q(i,M−i) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . (iii) For i∈{1, ..., M−1}, the sub-matrix Qi,i+1is specified by Qi,i+1= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ (1−p)βIi μi(1−p)βIi 2μi(1−p)βIi ...... (M−i−1)μi(1−p)βIi (M−i)μi ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . References Akar N, O˘guz NC, Sohraby K (2000) A novel computational method for solving finite QBD processes. Stoch Models 16:273–311 Allen LJS (2003) An introduction to stochastic processes with applications to biology. Pearson Education Inc, Upper Saddle River Almaraz E, Gómez-Corral A (2018) On SIR-models with Markov-modulated events: length of an outbreak, total size of the epidemic and number of secondary cases. Discrete Contin Dyn Syst Ser B 23:2153–2176 Artalejo JR, Lopez-Herrero MJ (2013) On the exact measure of disease spread in stochastic epidemic models. Bull Math Biol 75:1031–1050 Bhatta AK, Keyal U, Liu Y, Gellen E (2018) Vertical transmission of herpes simplex virus: an update. J Dtsch Dermatol Ges 16:685–692 Bichara D, Iggidr A, Sallet G (2014) Global analysis of multi-strains SIS, SIR and MSIR epidemic models. J Appl Math Comput 44:273–292 Busenberg S, Cooke K (1992) Vertically transmitted diseases. Springer, Berlin Chalub FACC, Gómez-Corral A, López-García M, Palacios-Rodríguez F (2023) A Markov chain model to investigate the spread of antibiotic-resistant bacteria in hospitals. arXiv:2305.14185 Chen Y, Zhao W (2020) Asymptotic behavior and threshold of a stochastic SIQS epidemic model with vertical transmission and Beddington-De Angelis incidence. Adv Differ Equ 2020:353 De Nitto Personè V, Grassi V (1996) Solution of finite QBD processes. J Appl Probab 33:1003–1010 D’Onofrio A (2005) On pulse vaccination strategy in the SIR epidemic model with vertical transmission. Appl Math Lett 18:729–732 Economou A, Gómez-Corral A, López-García M (2015) A stochastic SIS epidemic model with heterogeneous contacts. Phys A 421:78–97 Fernández-Peralta R, Gómez-Corral A (2021) A structured Markov chain model to investigate the effects of pre-exposure vaccines in tuberculosis control. J Theor Biol 509:110490 Gao S, Xie D, Chen L (2007) Pulse vaccination strategy in a delayed SIR epidemic model with vertical transmission. Discrete Contin Dyn Syst Ser B 7:77–86 Gaver DP, Jacobs PA, Latouche G (1984) Finite birth-and-death models in randomly changing environments. Adv Appl Probab 16:715–731 Gómez-Corral A, López-García M (2017) On SIR epidemic models with generally distributed infectious periods: number of secondary cases and probability of infection. Int J Biomath 10:1750024 Gómez-Corral A, López-García M (2018) Perturbation analysis in finite LD-QBD processes and applications to epidemic models. Numer Linear Algebra Appl 25:2160 Gómez-Corral A, López-García M, Lopez-Herrero MJ, Taipe D (2020) On first-passage times and sojourn times in finite QBD processes and their applications in epidemics. Mathematics 8:1718 Gómez-Corral A, Langwade J, López-García M, Molina-París C (2023) Sufficient conditions for regularity, positive recurrence and absorption in level-dependent QBD processes and related block-structured Markov chains. Math Methods Appl Sci 46:6756–6766 Gou Q, Wang W (2007) Global stability of two epidemic models. Discrete Contin Dyn Syst Ser B 8:333–345 123 291 Page 20 of 20 A. Gómez-Corral et al. He QM (2014) Fundamentals of matrix-analytic methods. Springer, New York Hu Z, Ma W, Ruan S (2012) Analysis of SIR epidemic models with nonlinear incidence rate and treatment. Math Biosci 238:12–20 Iglehart DL (1964) Multivariate competition processes. Ann Math Stat 35:350–361 Jaan A, Rajnik M (2021) TORCH Complex. In: StatPearls [Internet]. Treasure Island (FL): StatPearls Publishing. PMID: 32809363 Jeger MJ, van den Bosch Madden LV, Holt J (1998) A model for analysing plant-virus transmission characteristics and epidemic development. IMA J Math Appl Med Biol 15:1–18 Kang Y, Castillo-Chavez C (2014) Dynamics of SI models with both horizontal and vertical transmissions as well as Allee effects. Math Biosci 248:97–116 Kiouach D, Sabbar Y (2018) Stability and threshold of a stochastic SIRS epidemic model with vertical transmission and transfer from infectious to susceptible individuals. Discrete Dyn Nat Soc 2018:7570296 Kulkarni VG (2020) Modeling and analysis of stochastic systems, 3rd edn. CRC Press, Boca Raton Li MY, Smith HL, Wang L (2001) Global dynamics of an SEIR epidemic model with vertical transmission. SIAM J Appl Math 62:58–69 Li J, Blakeley D, Smith RJ (2011) The failure of R0. Comput Math Methods Med 2011:527610 Lu Z, Chi X, Chen L (2002) The effect of constant and pulse vaccination on SIR epidemic model with horizontal and vertical transmission. Math Comput Model 36:1039–1057 Ma JL, Earn DJD (2006) Generality of the final size formula for an epidemic of a newly invading infectious disease. Bull Math Biol 68:679–802 Marsico C, Kimberlin DW (2017) Congenital cytomegalovirus infection: advances and challenges in diagnosis, prevention and treatment. Ital J Pediatr 43:38 Naji RK, Hussien RM (2016) The dynamics of epidemic model with two types of infectious diseases and vertical transmission. J Appl Math 2016:4907964 Naresh R, Tripathi A, Omar S (2006) Modelling the spread of AIDS epidemic with vertical transmission. Appl Math Comput 178:262–272 Neu N, Duchon J, Zachariah P (2015) TORCH infections. Clin Perinatol 42:77–103 Reuter GEH (1961) Competition processes. In: Neyman J (ed) Proceedings of the 4th Berkeley symposium on mathematical statistics and probability, vol II: contributions to probability theory. University of California Press, Berkeley, pp 421–430 Roberts MG (2007) The pluses and minuses of R0. J R Soc Interface 4:946–961 Sullivan EM, Burgess MA, Forrest JM (1999) The epidemiology of rubella and congenital rubella in Australia, 1992 to 1997. Commun Disease Intell 23:209–214 van den Driessche P, Watmough J (2002) Reproduction number and sub-threshold endemic equilibria for compartmental models of disease transmission. Math Biosci 180:29–48 Wang X, Wang C, Wang K (2020) Global dynamics of a novel deterministic and stochastic SIR epidemic model with vertical transmission and media coverage. Adv Differ Equ 2020:685 Zhang XB, Chang S, Shi Q, Huo HF (2018) Qualitative study of a stochastic SIS epidemic model with vertical transmission. Phys A 505:805–817 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 123