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

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\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal{R}_{\mathrm{{exact}},0}$$\end{document}, which is seen to be the stochastic version of the well-known basic reproduction number. We characterize the probability distribution function of Rexact,0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal{R}_{\mathrm{{exact}},0}$$\end{document} 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.


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;Lu et al. 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.
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;Reuter 1961) on the positive orthant N 0 × N 0 , 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 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 R exact,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 R 0 is 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 R 0 > 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 R 0 overestimates the reproductive potential of a pathogen, especially when the involved community has a small size. An appropriate alternative to R 0 -in particular, for Markov chain models in epidemics-is the exact reproduction number R exact,0 , which is introduced in Artalejo and Lopez-Herrero (2013) by translating the definition of R 0 into a random variable. Among the most relevant features of R exact,0 , compared to R 0 , are that R exact,0 eliminates 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 R exact,0 requires knowing a priori the underlying parameters of the epidemic model. In computing the mass function of R exact,0 and its moments, the most intensive effort lies in the inversion of matrices and related memory requirements, whence the numerical solution of R exact,0 may experience instability with a large population size. These features and limitations show the adequacy of R exact,0 to quantify, in a more realistic way than R 0 , 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 number R H exact,0 of secondary infections by contact between susceptible individuals and a (predetermined) initially infectious individual, and the number R V exact,0 of 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 (R V exact,0 , R H 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.

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; see Fig. 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) = γ T min{i, C}, where γ T ≥ 0 denotes the treatment rate, 1 and C is 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β I or 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. The birth rates β S and β I are 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 = α I if the death of an individual does not depend on the compartment and, consequently, α I > α S when 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 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) , where q (i,s) is the rate of the sojourn time of process X in state (i, s); that is for (i, j) ∈ S, where δ a,b represents the Kronecker delta. Note that the state (0, 0) is absorbing, and X evolves as a birth-death process on the class {(0, s) : s ∈ {0, ..., M}} of infection-free states, with birth parameters {β S s : s ∈ {1, ..., M − 1}} and death parameters {α S s : s ∈ {1, ..., M}}.
For later use, we use a lexicographical ordering of states and decompose the state space This labeling of states allows to reformulate the process X as a finite QBD process (Gaver et al. 1984;Gómez-Corral et al. 2020 where the sub-matrices Q i,i are seen to be of dimension In particular, the non-vanishing entries of Q i,i−1 are related to the removal of an infectious individual due to its recovery or death; the diagonal entries of whereas the non-vanishing off-diagonal entries of Q i,i are linked to the birth and the death of a susceptible individual; and the non-vanishing entries of Q i,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 X will 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 X in 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 ∈ N 0 } 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 can be seen as sufficient for the first passage to the subset {(0, s) : s ∈ N 0 } of disease-free states to occur with certainty, and in a finite mean time, from any initial state (i, s) ∈ N 0 ×N 0 ; see, e.g., Gómez-Corral et al. (2023, Theorem 3) and Reuter (1961, Theorems 2-3). the vertical and the horizontal transmission of the disease, whence we decompose R exact,0 into two random contributions R V exact,0 and R H exact,0 according 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 R exact,0 into R V exact,0 and R H exact,0 is 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 R exact,0 -in terms of sample paths-or its expectation is greater than one, the values of R V exact,0 , and of R H 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 − p of infectious offspring (Table 1) is closely related to the disease, whence no action can be taken on the value of 1 − p in 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 * = D if the initially infectious individual dies before recovering; I * = N if it recovers due to natural causes; and I * = T if it recovers due to treatment. The aim is to characterize the joint probability law of (R V exact,0 , R H exact,0 ) by evaluating the conditional probabilities It is clear that R V exact,0 and R H exact,0 are inherently linked to events taking place during the time interval [0, τ ) and, consequently, the conditional probabilities P (1,s 0 ) (i * , r V , r H ), for status i * ∈ {D, N , T } and pairs (r V , r H ) ∈ N 0 × N 0 , satisfy the equalities These equalities are derived by observing that the left-hand side of (2)-(4) corresponds to the marginal mass function In evaluating the probability P (1,s 0 ) (i * , r V , r H ), for a predetermined status i * ∈ {D, N , T } and a fixed pair (r V , r H ) ∈ N 0 ×N 0 , we first define, for states (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 V infectious offspring and generates r H infectious contacts with susceptibles during the residual interval [t, τ ), and its status at time τ +0 is i * , given that (i, s) is the current state of X and 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 (0)

A first-step analysis of process X
We can derive a system of linear equations for the conditional probabilities P (i,s) (i * , r V , r H ), for (i, s) ∈ S\l(0), i * ∈ {D, N , T } and (r V , r H ) ∈ N 0 × N 0 , using the first-step analysis. To be concrete, it is found that 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) α I replaced by q −1 (i,s) γ N and q −1 (i,s) (1 − δ C,0 )γ T , respectively. For any status i * ∈ {D, N , T } and pairs (r V , r H ) ∈ N 0 × N 0 \{(0, 0)}, it is also seen that for states (i, j) ∈ S \ l(0). Equations (5) and (6) can be easily written in matrix form using column vectors p i (i * , r V , r H ), for i ∈ {1, ..., M}, with (1 + s)th entry P (i,s) (i * , r V , r H ) for integers s ∈ {0, ..., M − i}, and a suitable decomposition of sub-matrices Q i,i−1 and Q i,i+1 . Specifically, the rate sub-matrix Q i,i−1 is written in terms of where Q D i,i−1 , Q N i,i−1 , and Q T i,i−1 consist 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−1 is related to transitions associated with the death and recovery of the remaining i −1 infectious individuals. The rate sub-matrix Q i,i+1 is written in a similar manner as follows: where Q V i,i+1 and Q H i,i+1 record the infinitesimal transition rates associated, respectively, with the vertical and the horizontal secondary infections generated by the initially infectious individual before time τ , and Q * i,i+1 records the infinitesimal transition rates associated with secondary infections generated by the remaining i − 1 infectious individuals.
In matrix form, Eqs. (5) and (6) are given by for i * ∈ {D, N , T } and i ∈ {1, ..., M}, where 1 a is the column vector of order a of ones;

A general iterative scheme
It is important to note that, for a fixed pair (r V , r H ) ∈ N 0 × N 0 , we do not need solve (7) for i ∈ {1, ..., M}, in the unknown vectors p i (i * , r V , r H ), for the predetermined selec- Therefore, the application of block-Gaussian elimination to solve (8) results in where As a last remark, we note that, for a predetermined status i * ∈ {D, N , T } and a fixed pair (r V , r H ) ∈ N 0 × N 0 , Algorithm 1 evaluates the solution of the system of linear equations (5) and (6)  Step 1: j := 0; w := min{r V , r H }.

Numerical work and discussion
In this section, we present selected numerical results to explore the influence of the contact rate μ and the birth rate β I on the horizontal and vertical transmission of the disease during the early stages of an epidemic. We consider a population with maximum size M = 100, which consists initially of one infective and twenty susceptible individuals. In our experiments, three scenarios are defined in terms of the contact rate μ in such a way that less and more frequent infectious contacts occur in Scenarios 1 and 3, respectively, and are related to the values μ = 1.0 and μ = 5.0, and Scenario 2 corresponds to μ = 2.5. Birth and death rates are related to each other through the relationships β S = 2β I , α S = β S , and α I = 2α S , for values β I ∈ {1.25, 2.5, 3.75, 5.0, 10.0}. Therefore, the subpopulation of susceptible individuals is assumed to behave stable over time in terms of the underlying birth and death rates; in contrast, births among infectious individuals are expected to occur less frequently than births among susceptibles, and infectious individuals are more likely to die than susceptible ones. The value p = 0.2 is assumed, so that 80% of newborns from infectious parents are expected to be infective. It is assumed that an infectious individual recovers from natural causes, on In Figs. 4, 5 and 6, the marginal mass function of R exact,0 , the joint mass function of (R V exact,0 , R H exact,0 ) and the marginal mass functions of R V exact,0 and R H exact,0 , respectively, are plotted as a function of the contact rate μ and the birth rate β I . The specific intervals on the ox axis in Figs. 4 and 6, and on the axes ox and oy in Fig. 5 are selected to accumulate at least 99% of the probability of the corresponding probability distribution function. It is observed that, as intuition tells us, the higher the rate μ of contact between individuals, the more easily the pathogen spreads from the initially infectious individual to susceptible individuals at an early stage of the epidemic. This behavior is related to Fig. 4, where the marginal mass function of R exact,0 is seen to be more concentrated over the value r = 0 and values close to it for small values of the contact rate μ-as seen in Scenario 1, regardless of β I -, while it becomes more concentrated over higher values of r -for example, r = 4 in Scenario 3 for β I = 1.25-as μ increases. Unlike the behavior of the probability law of R exact,0 with respect to μ, increasing values of the birth rate β I do not lead to further spread of the pathogen in our numerical experiments; for instance, the mass function of R exact,0 in Scenario 1 becomes more concentrated on the value r = 0 and its surroundings as β I increases, contrary to what was expected. This apparently contradictory behavior can be explained by taking into account that births leading to infectious progeny from the initially infectious individual-that is, contributing to the value of R exact,0 -are recorded during the time interval [0, τ ) with expected length (α I + γ N + γ T ) −1 , where α I = 4β I . Thus, the length of this time interval is likely to decrease as β I increases and, consequently, fewer births of infectious individuals from the initially infectious individual contribute to R exact,0 , even though the birth rate β I is increasing. This makes the spread of the pathogen from the initially infectious individual to its progeny less likely with increasing values of β I . Figures 5 and 6 show, in more detail, how secondary infections due to horizontal transmission-in terms of R H exact,0 -and due to vertical transmission-in terms of R V exact,0contribute each one to the probability law of R exact,0 . 2 It is seen that the pathogen is more likely to be spread by person-to-person contact as fewer infectious progeny are given birth to by the initially infectious individual; i.e., on the set {R V exact,0 = r V } with r V ∈ {0, 1} in Fig.  5. Note that, in such a case, the value of R H exact,0 is expected to be higher as the number of person-to-person contacts increases. It is also remarkable that the number of person-to-person contacts has a notable influence on how the number of secondary infections due to vertical transmission can be predicted from the number of secondary infections due to horizontal transmission. This behavior can be observed by comparing Scenario 1 with Scenario 3 on In terms of the marginal probability laws of R H exact,0 and of R V exact,0 , Fig.6 shows that the marginal mass function of R H exact,0 is strongly influenced against the changes of μ and β I , while R V exact,0 is not very sensitive with respect to the variation of these parameters. In Fig. 6 The marginal mass functions of R H exact,0 (left) and of R V exact,0 (right) versus the contact rate μ and the birth rate β I , at an invasion time with s 0 = 20 particular, it is seen that the main contribution to R exact,0 in our experiments comes from the spread of the pathogen through person-to-person contacts.
To conclude, it is seen in Table 2 and Fig. 7 that the most likely value of R exact,0 -denoted by r max -behaves, as expected, as an increasing function of μ. Moreover, the corresponding conditional probability P(R exact,0 = r max |(I (0), S(0)) = (1, s 0 )) starts decreasing toward its minimum value and then increases with increasing values of μ, irrespectively of β I . On the contrary, the values of r max decrease with increasing values of β I in our numerical results, and the corresponding conditional probability does not have a common general behavior as a function of β I (Table 2).

Concluding remarks
In this paper, a Markov chain model to describe the spread of epidemics of SIS type through interactions between susceptible and infectious individuals, and generational transmission from colonized parents to their progeny has been presented. The analysis of the model requires the use of a bi-variate competition process X , which can be seen as the natural generalization of the uni-variate birth-death process arising in the classical SIS epidemic model without vertical transmission. Since bi-variate competition processes can be formulated as QBD processes, a broad range of mathematical techniques for block-structured Markov chains have great potential in the epidemiological setting to evaluate well-known descriptors; among these descriptors, we can mention the limiting distribution of the number of susceptible and infectious individuals as the time index tends to infinite, the random length of an outbreak, extreme values, and perturbation properties. In analyzing how the pathogen spreads in an early stage of the epidemic, the exact reproduction number R exact,0 has been expressed in terms of two random contributions allowing us to quantify how many secondary infections are related to the dynamics of horizontal and vertical transmission from an initially infectious individual. The joint probability law of the resulting random numbers R H exact,0 and R V exact,0 has been characterized in terms of finite systems of linear equations, which can be solved in an iterative manner by block-Gaussian elimination. It is clear that, since secondary cases are related to an exponentially distributed time τ , the probability law of R exact,0 , and those of R H exact,0 and R V exact,0 , are non-defective. This also holds in the case M = ∞, due to the regularity of the underlying QBD process with infinitely may states (Remark 1), and Eq. (8) is seen to remain valid for integers i ∈ N.
A major problem relates to the understanding of mechanisms which are responsible for the dissemination of a pathogen at the person-to-person and generational level. The specialized versions R H exact,0 and R V exact,0 are useful to understand how, specially in the case of populations with stable size over time, the person-to-person behavior can influence the dynamics at an early stage of the epidemic more strongly than the generational transmission from a initially colonized parent to its progeny. Therefore, the decomposition of the exact reproduction number into R H exact,0 and R V exact,0 can be considered as a promising tool for the study of how the mechanisms of vertical transmission are important in modeling the epidemic. Under the particular circumstances that the marginal mass function of R V exact,0 is seen to be essentially concentrated over the value r V = 0, the bi-variate process X could be appropriately replaced by a simple uni-variate birth-death process and, consequently, the corresponding analysis of related descriptors might be analytically simplified.
Other Markov chain models for epidemics of the SIS type with vertical transmission may be appropriately studied by adapting our methodology. By way of example, we con-  (Fig. 8), and we exemplify how the QBD formulation is the key in the derivation of the joint distribution of (R V exact,0 , R H 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 ): ..., M}}, and q (i,s,v), (i,s,v) for (i, s, v) ∈ S . In these expressions, α V is 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. Probabilitiesp andp , withp,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 probabilityp , withp ∈ [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 , the q-matrix of X is seen to have the structured form of Q in Eq. (1) using the labeling of states , ..., M}. This means that, for process X , sub-matrices Q i,i−1 and Q i,i+1 in (1) record infinitesimal rates associated with the recovery and the death of infectious individuals, and new infections, respectively, and diagonal entries of Q i,i are specified by −q (i,s,v) , for (i, s, v) with s + v ∈ {0, ..., M − i}. Non-vanishing off-diagonal entries of Q i,i correspond to the birth and the death of susceptible and vaccinated individuals, and the therapeutic use of a vaccine on susceptible individuals. The column vectors p i (i * , r V , r H ) for process X , for integers i ∈ {1, ..., M} and i * ∈ {D, N , T }, are defined as p i (i * , r V , r H ) in Section 3.1 with the probabilities P (i,s) (i * , r V , r H ) replaced by the conditional probabilities P( (9) and (10) are then found to determine the joint probability law of (R V exact,0 , R H exact,0 ) by adapting our approach in Sections 3.1 and 3.2 to process X ; in particular, matrices Q i,i−1 (respectively, Q i,i+1 ) for X 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 {R V exact,0 ∈ {0, ..., r V }} and {R H exact,0 ∈ {0, ..., r H }}, for certain relatively small values r V and r H , are expected to occur more likely in process X than in process X and, on at least first observation, the tail of the distribution of R V exact,0 (respectively, R H exact,0 ) for process X is expected to be heavier than that of R V exact,0 (respectively, R H 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 ofp,p andp -compared to a therapeutic intervention-defined by a strictly positive value of λ andp =p =p = 0especially to derive certain thresholds for the underlying probabilitiesp,p andp , and the 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 X can be defined by replacing the exponentially distributed recovery times with mean γ −1 N by 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 = (y 1 , ..., y m )-with m j=1 y j = i-recording the number y j of 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 (R V exact,0 , R H exact,0 ) will become rapidly prone to numerical 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).