MD A J Van Sorge¹

student PhD J U M Termote²

PhD M J De Vries

MD F N Boonstra

PhD C Stellingwerf

MD N E Schalij-Delfos¹

MD, Student-PhD A J Van Sorge
email: j.van_sorge@lumc.nl

The incidence of visual impairment due to ROP and concomitant disabilities in the Netherlands

Keywords:

Purpose: Determine the incidence of visual impairment (VI) caused by Retinopathy of Prematurity (ROP) and the concomitant disabilities in preterm neonates born between 2000 and 2009 in The Netherlands.

Methods: Data were retrieved from the Dutch institutes for the visually impaired. They were compared to similar Dutch studies conducted between 1975-1987, 1986-1994 and 1994-2000. Results: Records of 42 infants with VI due to ROP were included. A gradual decrease of gestational age and birth weight but an increase of duration of artificial ventilation, supplemental oxygen administration, bronchopulmonary dysplasia, developmental delay and behavioral abnormalities was found. Compared to the previous study (1994-2000) significantly less children were visually impaired due to ROP (1.84 per 100.000 live births/year vs 3.93 per 100.000 live births/year, p=0.000), the incidence of complete blindness decreased from 27.5% to 7.1% (p<0.05) and more children were treated (66.7% vs 56.9%, ns).The incidence of concomitant disabilities was again high and did not differ a lot from the previous study.

Conclusion:

Retrospective study showing a significant decrease in VI due to ROP in the Netherlands.

Changes in neonatal care did not result in a decrease of concomitant disabilities. More children were treated for ROP, but still 33% were not.

Introduction

Retinopathy of prematurity (ROP) is still one of the most important causes of partial sight or blindness in premature born infants. In developed countries, ROP accounts for 5.5 to 20% of childhood blindness. [START_REF] Schiariti | Severe retinopathy of prematurity and visual outcomes in British Columbia: a 10-year analysis[END_REF][START_REF] Gilbert | Retinopathy of prematurity: a global perspective of the epidemics, population of babies at risk and implications for control[END_REF][START_REF] Khan | Changing pattern of childhood blindness[END_REF][START_REF] Rahi | Severe visual impairment and blindness in children in the UK[END_REF] The pathogenesis of ROP is multifactorial, but the most significant risk factors for ROP are low gestational age and low birth weight. [START_REF] Slidsborg | Treatment for retinopathy of prematurity in Denmark in a ten-year period (1996 2005): is the incidence increasing?[END_REF][START_REF] Termote | Can screening for retinopathy of prematurity be reduced?[END_REF][START_REF] Dhaliwal | Retinopathy of prematurity in smallfor-gestational age infants compared to appropriate-for-gestational age infants[END_REF] Changes in neonatal care over the past two decades have increased survival of very small and very preterm born infants. This shift was due to introduction of surfactant, increased use of antenatal steroids, resuscitation of infants with lower gestational age, improved nutrition and prophylactic use of indomethacin. Due to this increased survival, a larger number of very preterm infants is at risk not only for visual impairment (VI) caused by ROP but also for non-visual handicaps. [START_REF] Hack | Survival and neurodevelopmental outcomes of preterm infants[END_REF][START_REF] Msall | Measuring functional outcomes after prematurity: developmental impact of very low birth weight and extremely low birth weight status on childhood disability[END_REF][START_REF] Colvin | Neurodevelopmental outcomes after preterm birth[END_REF] Apart from changes in neonatal care, there have been changes in management of ROP. To reduce the number of adverse outcomes, treatment in the pre-threshold phase was introduced, resulting in an increased number of infants treated for ROP. 5;11;12 From 1952 onwards, five surveys were carried out in the Netherlands to calculate the number of

Patients and Methods

Due to privacy regulations and lack of an obligatory registry for children with VI in our country, the access to patient data is limited. Therefore, data retrieval was only possible by collecting the records of all infants with a registered diagnosis of ROP, who are known at one of the Dutch institutes for the partially sighted and blind. This information was obtained by one investigator (AvS).

For comparison the same neonatal and ophthalmological data as in the previous surveys were collected: gestational age (GA), birth weight (BW), visual acuity, treatment for ROP, sex, multiple birth, incidence of bronchopulmonary dysplasia (BPD) and duration of supplementary oxygen administration, artificial ventilation (AV) and admission to the Neonatal Intensive Care Unit (NICU).

With regard to concomitant disabilities the same definitions and classifications as described by Termote were used. [START_REF] Termote | The incidence of visually impaired children with retinopathy of prematurity and their concomitant disabilities[END_REF] We considered children to be multiple disabled when they had VI caused by ROP plus one or more of the concomitant disabilities that were listed, excluding BPD as pulmonary function usually improves through the years.

To define VI, the recommendations of the International Association for Prevention of Blindness (IAPB)

were used (WHO 1984). Due to limitations in data collection, information on stages of ROP was unreliable and therefore again excluded. As the number of infants with cerebral visual impairment (CVI) increased over the last decade and became the most important cause of VI in prematurely born infants [START_REF] Good | Cortical visual impairment: new directions[END_REF] , data on CVI were collected for this survey. They could however not be compared to the previous studies.

Dutch birth rate figures were obtained from the Central Bureau of Statistics (CBR). Survival rates of premature infants according to gestational age were collected from the Dutch Neonatal Registration.

The study was approved by the Medical Ethical Committees of all participating institutes for the partially sighted and blind and children were included when informed consent was obtained.

Data collected for this study (2000-2009) were compared to data from the previous study (1994-2000).

Besides, data of the four surveys conducted from 1975 onwards were compared to evaluate which trends in treatment and outcome could be seen.

Statistical analysis

Results

Records of 43 children with VI caused by ROP were found. Informed consent was obtained from 42 parents. No additional infants born between 1994 to 2000 (survey 5) registered with the diagnosis of ROP were found during survey 6.

General and neonatal data of infants diagnosed with ROP in four consecutive periods in The Netherlands are presented in table 1. Compared to the previous study (1994-2000) no significant differences were found. However, a gradual decrease of mean GA and BW over the years was seen.

For supplemental oxygen and AV a significant increase was found between survey 3 and 4, but there were no significant changes since the mid-eighties. Not all infants from survey 3 and 4 were admitted to a NICU, so admission days to a NICU were only analyzed for the fifth and sixth period, for which there was no significant difference. A decline can be seen starting in the mid-nineties. The incidence of VI due to ROP decreased in the period 1994-2000 but failed to reach significance when compared to period 1986-1994 (p=0.07).

Comparing our study to survey 4 and 5 a significant decrease of the incidence of VI caused by ROP was found (p=0.000 for both). More specific data on visual acuity are presented in table 3. In the current study there were 32 of 42 infants (76.2%) with VI, compared to 46 of 51 children (90.2%) in survey 5, representing a non significant decrease. In more detail, our survey shows a significant decrease of completely blind children due to ROP compared to the previous period (7.1% vs 27.5%; p=0.012) as well as a significant increase of children who did not fulfill the WHO criteria for VI at the moment of inclusion in this study (23.8% vs 9.8%; p=0.043). Eight children (19%) were diagnosed with CVI: 2 not partially sighted or blind, 4 partially sighted, 1 socially and 1 practically blind. The number of children who were treated for ROP significantly increased since period three, but there was no significant increase between period 5 and 6.

Infants were treated with laser (n=16), cryo (n=1), laser + cryo (n=2), cerclage (n=1), cerclage + laser (n=1), cerclage + cryo (n=1), vitrectomy + lensectomy (n=2), vitrectomy + laser (n=2), vitrectomy + laser + lensectomy (n=1) or bevacizumab (n=1). Most common therapy used was laser (79%).

Fourteen infants (33.3%) of the present study had not been treated against 22 (43.1%) in the previous study. Four of the untreated children in this study had a final VA > 0.3 at least in one eye.

Comparing the incidence of concomitant disabilities in infants with ROP, no significant changes were found except for epilepsy (p=0.007) (table 4).

Table 4: Concomitant disabilities in infants with visual impairment caused by ROP in four consecutive periods in The Netherlands (%).

¹ period 1 vs 2 p=0.009, 1 vs 3 p=0.000, 3 vs 4 p= 0.016 ²1 vs 2 p= 0.028, 1 vs 3 p= 0.000, 1 vs 4 p= 0.000, 2 vs 3 p= 0.000, 2vs 4 p= 0.027, 3 vs 4 ns ³1 vs 3 p= 0.049 and 3 vs 4 p=0.007 4 2 v s 4 p = 0 . 0 2 4 [START_REF] Slidsborg | Treatment for retinopathy of prematurity in Denmark in a ten-year period (1996 2005): is the incidence increasing?[END_REF] 1 vs 4 p=0.003 [START_REF] Termote | Can screening for retinopathy of prematurity be reduced?[END_REF] 1 vs 2 p=0.013 7 1 vs 2 p=0.015 , 1 vs 3 p=0.0011 vs 4 p=0.001

In the previous study a significant increase in behavioral abnormalities was seen, but this tendency did not extend to the current study. A gradual increase of with BPD, developmental delay and hearing deficit can be observed, but differences did not reach significance. The number of infants with neurological handicaps and a multiple disability increased significantly between survey 3 and 4 but was quite stable afterwards.

When the population was split up for GA, the incidence of VI caused by ROP based on the estimated number of survivors as well as the incidence of concomitant disabilities in period 5 and 6, was the highest in children under 30 weeks and concerning VI, increased as GA decreased (table 5). The five infants with unilateral blindness (n=1) or VA > 0.3 (n=4) were also included in this analysis, because this was also done in the previous studies. The percentage of children with concomitant disabilities was again high in this survey, but did not differ a lot from survey number 5. In survey 5 there were 35 children with concomitant disabilities (68.6%) versus 31 (73.8%) in our survey.

Discussion

We performed a retrospective study on the incidence of VI due to ROP in preterm infants born between 2000 and 2009 in the Netherlands. The incidence presented in this study must be considered as a minimum incidence because there is no obligatory national registry for VI, children may be missed and there is always a possibility that children born in the study period will become visually impaired at a later age. However the same applied for earlier studies and no additional infants of the previous periods were found during our current search. Data were compared with three earlier Dutch studies starting in 1975 (survey 3, 4 and 5), thus providing an overview of more than thirty years.

Over the years there has been an increase of survival of smaller and more immature infants due to changes in neonatal care. [START_REF] Valcamonico | Mid-and long-term outcome of extremely low birth weight (ELBW) infants: an analysis of prognostic factors[END_REF] From 1975 till now mean GA and BW of infants with VI due to ROP in the Netherlands gradually decreased. Although we found no significant changes in neonatal data in infants with VI in the present study compared with the former one, we found a decrease in VI due to steroids and surfactant and a decrease in duration of AV, supplemental O 2 administration and NICU admission as well as a decrease in incidence of ROP. Most of these changes were not reproducible in our study. However, we included only those infants that were referred to one of the Dutch institutes for the visually impaired and were supposedly the most critically ill neonates. We speculate that more subtle improvements in neonatal care are responsible for the decrease of the incidence of VI caused by ROP and of completely blind children. Furthermore there was an increase in infants treated for ROP (from 56.9 to 66.7%). The Early Treatment criteria (ETROP) were introduced in the early phase of this survey resulting in a larger number to treat and an expected reduction of adverse outcomes. 3;11;12 Besides, the former study resulted in a campaign for more awareness of necessity to treat timely and for potential visual problems on the long term. Earlier referral to the institutes together with changed treatment criteria could be an additional explanation for the significantly higher number of infants with a VA > 0.3 found in the current study and the significant decrease in VI due to ROP.

Although more infants received treatment in the current study, still 33 % did not of which 10 had a VA < 0.3. Haines et al [START_REF] Haines | UK population based study of severe retinopathy of prematurity: screening, treatment, and outcome[END_REF] also distinctly described that, even after implementation of screening guidelines, some children developed VI due to non-adherence to the screening protocol.

Studies from the UK on VI due to ROP as a proportion of childhood VI, covering the same period 4;23 show a tendency that compares favorably with our surveys. A study of Schiariti et al This overview over more than 30 years shows a decrease in GA and BW and a gradual increase in AV and supplemental oxygen administration with a turning point in the early nineties when novel treatments such as surfactant treatment, high frequency oscillation (HFO), postnatal steroids and inhaled NO were introduced. The incidence of VI caused by ROP decreased over the past thirty years with the most obvious change seen in the last period. Additional changes in neonatal care such as antenatal steroids and antenatal antibiotics together with earlier treatment for ROP could be an explanation. However, these changes did not result in a decrease of concomitant as well as multiple disabilities.

Although more children were treated, still 33% were not. Therefore, it can not be emphasized enough that timely screening, follow up and treatment as well as close collaboration of all professionals involved in the care of these at risk neonates should be our aim for the future to further decrease the number of children with VI caused by ROP. As long as GA and BW of survivors decline, it seems a more difficult task to diminish the number of multiple disabled children.

  visually disabled children due to ROP per 100.000 live births per year: 3.24/100.000 from 1952-1964 [Schappert-Kimmijser; survey 1], 2.89/100.000 from 1961-1973 [vd Pol; survey 2], 4.22/100.000 from 1975-1987 [Cats; survey 3], 5.49/100.000 from 1986-1994 [Schalij-Delfos; survey 4] and 3.85/100.000 from 1994-2000 [Termote; survey 5]. The last 3 surveys also included an inventory of concomitant disabilities, showing a marked increase of developmental and behavioral abnormalities. To monitor trends and to determine the incidence of VI caused by ROP as well as the incidence of associated disabilities in these high-risk, vision-impaired children born between 2000 and 2009, a retrospective study was performed [survey 6]. Data of the last 3 surveys were used for comparison, resulting in an overview of more than three decades.

  ROP. It seems reasonable to assume that improvement in neonatal care might contribute to a better visual outcome of these preterm neonates. An important limitation is that we could not have more insight in detailed information about the neonatal period. From other Dutch studies covering approximately the same period we can learn more about general changes that have taken place. Two Dutch cohort studies from Hoogerwerf[START_REF] Hoogerwerf | Incidence of Retinopathy of Prematurity over the Last Decade in the Central Netherlands[END_REF] and Groenendaal[START_REF] Groenendaal | Complications affecting preterm neonates from 1991 to 2006: what have we gained?[END_REF] demonstrated an improvement in outcome. Groenendaal et al analyzed two cohorts of inborn preterm neonates (1997-2001) and (2002-2006) with a GA of 25-29.9 wks and found an improved survival but no increase in the number of infants treated for severe ROP. The study from Hoogerwerf et al focused on changes in the incidence and risk factors responsible for the development of ROP and compared two cohorts of prematures born in the time period of 1991-1995 and 2001-2005. They found an increase in the use of antenatal

  Clinical data of the four surveys as presented in table 1 were evaluated using the

	independent samples t-test. The Poison regression was used to calculate figures on VI caused by
	ROP in relation to Dutch birthrates and the chi-square test for all other data. Differences with a p-value
	of <0.05 were considered significant. As all children with VI due to ROP were studied in the
	subsequent periods, trends could be determined even though the separate values failed to reach
	significance.

Table 1 : Neonatal data of infants with visual impairment caused by ROP in four consecutive periods in The Netherlands.

 1

	Period / survey #		1975-1987 / 3	1986-1994 / 4	1994-2000 / 5	2000-2009 / 6
	No. Infants		76	87	51	42
	Male (%)		47	58	59	74
	Mean gestational age (wks)¹	28.7 ± 2.7	27.5 ± 2.5	27.7 ± 2.4	27.4 ± 2.0
	Mean birth weight (gr)²		1128 ± 331	1071 ± 385	942 ± 306	912 ± 385
	Multiple birth (%)		18.9	20.7	31.4	31.0
	Supplemental O₂ administration	43.2 ± 40.7	78.6 ± 88.8	83.2 ±81.6	81.8 ±78.4
	(days)³
	Artificial ventilation (days)	4	13.3 ± 22.0	27.6 ± 27.7	24.6 ±16.6	25.3±17.5
	Admission to NICU (days)	5	-	-	92.1 ± 58.2	97.2 ± 71.1

¹ column 1 vs column 2 p= 0,004, 1 vs 4 p=0.004, 2 vs 3, 3 vs 4 and 2 vs 4 ns ² column 1 vs 2 ns, 1 vs 3 p=0.001, 1vs 4 p=0.003, 2 vs 3 p=0.032, 3 vs 4 ns, 2 vs 4 p=0.03 ³ column 1 vs 2 p=0.003, 1 vs 3 p=0.005, 1 vs 4 p =0.004, 2 vs 3, 3 vs 4 and 2 vs 4 ns 4 column 1 vs 2 p=0.0009, 1 vs 3 p=0.004, 1 vs 4 p=0.004, 2 vs 3, 3 vs 4 and 2 vs 4 ns 5 column 3 vs 4 ns Although no significant differences were found in clinical data of children from the present study as compared to the previous study, the proportion of infants with VI caused by ROP in relation to Dutch birthrates decreased significantly (table

2

).

Table 2 : Visual impairment caused by ROP in relation to Dutch birth rates.

 2

*Central Bureau of Statistics for the Netherlands

The average incidence of VI due to ROP is presented column 3 is the result of the fraction of the absolute number of individuals with ROP sequelae (column 1) and the absolute number of individuals born during the study period (column 2). The average incidence of VI due to ROP shows a decrease since the early nineties.

Table 3 : Visual acuity in children with ROP in four consecutive periods

 3

	Period		1975 -1987 survey 3	1986 -1994 survey 4	1994 -2000 survey 5	2000 -2009 survey 6
	No. Infants	76	87	51	42
	Visual acuity unspecified	5.1%
	Not partially sighted or blind
	(VA>0.3)	+ ¹	2%	10.3%	9.8%	23.8%
	Partially sighted (VA 0.1 -0.3)	34.3%	31.0%	25.5%	38.1%
	Socially blind (VA < 0.1 -≥ 1/60)	12.1%	10.3%	11.8%	14.3%
	Practically blind (VA < 1/60 -LP)	8.1%	21.8%	25.5%	16.7%
	Completely blind (VA = 0)²	38.4%	26.4%	27.5%	7.1%
	Treatment of acute ROP³	24.%%	43.9%	56.9%	66.7%
	LP= light perception
	V A = v i s u a l a c u i t y
	¹ 1 vs4 p= 0.000, 2 vs 4 p= 0.043
	² 1 vs 4 p=0.000, 2 vs 4 p=0.010 and 3 vs4 p=0.012
	³ 1 vs 2 p=0.013, 1 vs 3 p=0.000, 1 vs 4 p=0.000, 2 vs 4 p=0.014

+ = infants with unilateral blindness and infants whose vision is expected to deteriorate in the future

Table 5 : Visual impairment caused by ROP in relation to estimated number of survivors and concomitant disabilities.

 5

	GA (wks)	% VI caused by ROP (n) related to estimated survivors	% Concomitant disabilities (n)
		'75 -'87	'86 -'94	'94 -'00	'00 -'09	'94 -'00	'00 -'09
		survey 3	survey 4	survey 5	survey 6	survey 5	survey 6
	< 24
	24 -25	10.4 (11)	12.5 (11)	3.1 (9)	2.4(10)	67 (6)	70 (7)
	26 -27	3.3 (40)	4.3 (36)	1.3 (23)	0.9(19)	65 (15)	79 (15)
	28 -29	0.8 (26)	0.4 (10)	0.3 (10)	0.2(9)	80 (8)	78 (7)
	30 -31	0.2 (13)	0.2 (8)	0.1 (6)	0.03(2)	66 (4)	50 (1)
	>32 -<37	0.01 (12)	0.01 (7)	<0.001 (2)	<0.002 ( 2)	100 (2)	50 (1)

* Numbers of estimated survivors were derived from The Netherlands Perinatal Registry

  Slidsborg et al[START_REF] Slidsborg | Treatment for retinopathy of prematurity in Denmark in a ten-year period (1996 2005): is the incidence increasing?[END_REF] found a more than twofold increase in the incidence of treated ROP cases born in Denmark from 2001 to 2005, compared with 1996-2000. This increase was most pronounced for the smallest infants. The number of infants with VI or blindness caused by ROP was stable. These studies also attribute the decline of VI due to ROP to a more complete ophthalmic surveillance, earlier treatment and improvement in neonatal care.Compared to the earlier Dutch surveys, again most children with VI caused by ROP were surviving preterm neonates < 30 wks of gestation. This group of infants also showed the highest incidence of concomitant disabilities. Although a substantial number of infants had 1 or more concomitant disabilities (73.8%), no significant changes were found compared to the study of Termote.Inexplicably, no children with epilepsy were registered in the present study. BPD, developmental delay and hearing deficit gradually increased whereas behavioral abnormalities, multiple disabilities and neurological handicaps decreased minimally. A Dutch study of van Baar et al[START_REF] Van Baar | Very preterm birth is associated with disabilities in multiple developmental domains[END_REF] showed a clear association between preterm birth and multiple disabilities at 5.5 years of age in 157 children born before 30 weeks GA. One or more disabilities were found in 75%. Several other studies conclude that the survival rate for extremely preterm born infants has improved over the past decade, but that the overall prevalence of neurodisability after preterm birth has not fallen thus supporting the findings in our study.

[START_REF] Schiariti | Severe retinopathy of prematurity and visual outcomes in British Columbia: a 10-year analysis[END_REF]

, demonstrated an increase of ROP as well as severe ROP without changes in the rate of VI when comparing neonates born in

1992-1996 and 1997-2001. 8;9;10


J R Artalejo
email: jesusartalejo@mat.ucm.es

M J Lopez-Herrero

Quasi-stationary and ratio of expectations distributions: a comparative study

Introduction

There exists a vast literature that uses Markov chains to model a biological population. A good reference to find the underlying mathematical theory is the textbook by [START_REF] Allen | An Introduction to Stochastic Processes with Applications to Biology[END_REF] and its references. As examples of recent applications of continuous-time Markov chains in population biology, we mention some papers published in this journal [START_REF] Casagrandi | The intermediate dispersal principle in spatially explicit metapopulations[END_REF][START_REF] Fouchet | Impact of myxomatosis in relation to local persistence in wild rabbit populations: The role of waning immunity and the reproductive period[END_REF][START_REF] Stirk | Stochastic niche structure and diversity maintenance in the T cell repertoire[END_REF][START_REF] Xu | Stochastic model of an influenza epidemic with drug resistance[END_REF].

Stochastic biological models based on Markov chains typically consist of a set S T of transient states where the process evolves until it escapes to a set 1 of absorbing states S A . In many applications S A reduces to the sole state 0. Due to the reducible character of the Markov chain, the stationary distribution becomes degenerate in the sense that it only assigns probability mass to the absorbing states. This fact motivates the following fundamental problem:

There exist a need for probabilistic measures of the system behavior before absorption. Analogues of the stationary distribution of the irreducible case should be considered and compared.

The above description of the problem under study reflects the spirit of the seminal work by Darroch andSeneta (1965, 1967) who suggested several possibilities including quasi-stationary distributions and ratio of means.

The term quasi-stationarity refers to the distribution of the Markov chain by conditioning on the event that absorption has not occurred yet. Let us assume that the Markov chain starts with the initial distribution on the transient states {u i ; i ∈ S T }. Then, the distribution {u i ; i ∈ S T } is called a quasistationary distribution when the distribution of the Markov chain, conditional on it has not yet absorbed at time t, remains equal to {u i ; i ∈ S T } over time. On the other hand, the quasi-stationary distribution can also be presented as a limiting conditional probability [START_REF] Darroch | On quasi-stationary distributions in absorbing continuous-time finite Markov chains[END_REF][START_REF] Pollett | Quasi-stationarity in populations that are subject to largescale mortality or emigration[END_REF]. This interpretation motivates the significant role of the quasi-stationarity as a measure of the time spent by the process in the transient states when the time to absorption is sufficiently long. For a formal definition of quasi-stationarity, we refer the reader to Section 2.

If the transient set S T is finite and irreducible, then the existence of the quasi-stationary distribution is guaranteed. However, even in the case of birth and death processes, it is typically impossible to evaluate explicitly the quasistationary distribution. One is therefore lead to consider iterative methods, approximations and asymptotic results for the quasi-stationary distribution. Many satisfactory approximate results have been obtained, see Nåsell's papers (1996Nåsell's papers ( , 1999aNåsell's papers ( , 2001) and their references. There are important models in population biology where the transient set S T is reducible. For a comprehensive treatment of this case, see [START_REF] Van Doorn | Survival in a quasi-death process[END_REF]. Roughly speaking, one could say that the quasi-stationary distribution, if it exists, gives positive probability to a certain subset S T + of S T but the quasi-stationary probability of a state i ∈ S T -S T + is zero. This is a striking fact which merits some extra attention. An example of this situation is given by the SIR epidemic model, see Subsection 3.2.1. If S T is infinite, the analysis is essentially more complicated. Conditions for the existence of the quasi-stationary distribution are mathematically involved. An attempt to summarize conditions given in terms of the infinitesimal transition rates can be found in [START_REF] Pollett | Quasi-stationarity in populations that are subject to largescale mortality or emigration[END_REF]. The book by [START_REF] Seneta | Non-Negative Matrices and Markov Chains[END_REF] is a basic reference for the general mathematical theory.

The above discussion on the concept of quasi-stationarity, which will be extended in Section 2, gave the authors an initial motivation to retrieve the RE-distribution (ratio of expectations distribution) as another natural measure of the behavior of absorbing Markov chains before absorption occurs. Given the initial state, the RE-distribution will be defined as a ratio between the time that the process spends at each transient state and the expected time to absorption. In the sequel, we will show how the RE-distribution has a meaningful interpretation. Moreover, it exists and can be evaluated with a great generality. The idea of presenting a ratio of means as an analogue of the stationary distribution of an irreducible Markov chain was suggested in the early work by [START_REF] Ewens | The diffusion equation and a pseudo-distribution in genetics[END_REF] and Darroch andSeneta (1965, 1967).

One of the fundamental reasons to study biological models is to provide tools for getting insight of the underlying biological reality. Different models and approaches can be considered, but the final objective is always to capture the essential features of a biological system. It would be difficult and subjective to decide which approach is preferable. In fact, stochastic models are based on the concept of probability, whose basic roots are indeed subject to different interpretations (see e.g. the different methods for quantifying uncertainty or the different concepts of convergence of random variables). The existence of several methods and approaches should be interpreted as an advantage rather than as a problem. This open spirit inspires our comparative study in the paper. It is clear that the quasi-stationarity is a well established concept and there is no doubt about its usefulness. Our aim is to rescue the ratio of expectations approach and show that it is worthy of some extra attention. The comparison between the quasi-stationary distribution and the RE-distribution is justified only if the convergence to quasi-stationarity is relatively fast. If this is the case, the comparative analysis done in the sequel may lead to the conclusion that both approaches give similar or different information. As a result, one may see pros and cons. However, it is so far of our intention to suggest the possible superiority of one approach over the other. The existence of different alternatives to measure the system behavior before absorption enriches our knowledge and provides complementary insight of the underlying biological problem.

The rest of the paper is organized as follows. In Section 2, the RE-distribution is introduced in the context of a general absorbing continuous-time Markov chain. Some emphasis is done on a general comparison with respect to the quasi-stationary distribution. Section 3 explores applications to stochastic population models. More specifically, we first consider the irreducible unidimensional case and deal with a class of birth and death processes including the stochastic SIS epidemic model (Subsection 3.1.1) and a linear growth model (Subsection 3.1.2). The SIR epidemic model (Subsection 3.2.1) provides a nice example to deal with the reducible case. As an example of a bidimensional tractable case where S T is irreducible, we consider a finite SIR model with demography (Subsection 3.2.2). These models provide appropriate scenarios to illustrate the tractability and interest of the RE-distribution as a measure of the system behavior before absorption. Moreover, a comparison with results based on the quasi-stationary distribution is performed. Our comparative study includes some theoretical justifications and discussion on the related computational issues. Finally, in the discussion section (Section 4), the main features of the quasi-stationary distribution and the RE-distribution are recapitulated.

The discrete-time case is also briefly outlined.

Comparing quasi-stationarity and the RE-distribution

We shall employ an absorbing continuous-time Markov chain {X(t); t ≥ 0} with countable state space S to model a biological population. Suppose that S A denotes the set of absorbing states and S T is the set of transient states, so

S = S A ∪ S T . Let Q = [q ij ]
be the infinitesimal generator of the Markov chain. Then, q i = -q ii represents the rate of the exponential sojourn time at state i. For states i ∈ S A , we define q i = 0.

The quasi-stationary distributions have been widely used to measure the distribution of the process given that the absorption has no yet taken place. Thus, the starting point is the conditional probability

u i (t) = P {X(t) = i | T > t}, i ∈ S T ,
where T = sup{t ≥ 0: X(t) ∈ S T } denotes the absorption time.

Now suppose the Markov chain starts with the initial distribution a i = P {X(0) = i}, for i ∈ S T . If there exists a starting distribution a i = u i such that

P {X(t) = i | T > t} = u i , i ∈ S T , (2.1) for all t ≥ 0, then {u i ; i ∈ S T } is called a quasi-stationary distribution.
From definition (2.1), it is clear that the quasi-stationary distribution is independent of t if the Markov chain starts with initial distribution {u i ; i ∈ S T }.

Let us first assume that S T is finite and irreducible. Then, it is wellknown [START_REF] Darroch | On quasi-stationary distributions in absorbing continuous-time finite Markov chains[END_REF] that the quasi-stationary distribution is a left eigenvector corresponding to the eigenvalue with maximal real part of the sub-generator Q S T associated to the transient states. We denote such maximal eigenvalue as -α. It is well-known that α is real and positive. This result provides a first option for numerical computations.

There also exists a limiting interpretation of the quasi-stationary distribution which states that

lim t→∞ P {X(t) = i | T > t} = u i , i ∈ S T , (2.2) independently of the initial distribution {a i ; i ∈ S T }.
The limiting result (2.2) supports the accepted idea that the quasi-stationary distribution provides a good measure of the system dynamics before absorption, but restricting only to those realizations in which the time to absorption is sufficiently large. Another doubly-limiting approach to quasi-stationarity has also been considered in the literature, see for instance [START_REF] Darroch | On quasi-stationary distributions in absorbing continuous-time finite Markov chains[END_REF] and [START_REF] Pollett | Quasi-stationarity in populations that are subject to largescale mortality or emigration[END_REF].

There is no doubt about the useful information provided by the quasistationary distribution in the case where the Markov chain is known to be in the set S T after a very long time. However, there exists a second scenario in which the quasi-stationary distribution is also helpful to describe the behavior of the system before absorption. This is the case where the process reaches the quasi-stationary regime after a relatively short time and, after a substantially longer period, absorption will occur.

At this point, a rule to identify the second scenario in practice is needed. A natural solution is given in terms of a comparison between the rate of convergence to absorption and the rate of convergence to the quasi-stationary distribution. The former is given by α, while the latter is αα where -α is the real part of the eigenvalue with the next smallest real part, see formulas (1.1) and (3.3) in [START_REF] Darroch | On quasi-stationary distributions in absorbing continuous-time finite Markov chains[END_REF]. As a result, the quasi-stationary distribution should be used to describe the behavior before absorption only if R u = 2α/α < 1. Then, the comparison between the quasi-stationary distribution and the RE-distribution becomes meaningful.

It is well reported in the literature [START_REF] Nåsell | Extinction and quasi-stationarity in the Verhulst logistic model[END_REF]) that an analytical solution to the quasi-stationary probabilities only exists in a few special cases. Among them, we mention some birth and death processes where the birth and death rates are either constants or linear functions of the state i ∈ S T . As a result, approximations and numerical methods of computation are required.

If we focus on finite birth and death processes with absorbing state 0, then the difficulties to obtain an analytical solution are related to the non-linear form of the system of equations governing the quasi-stationary probabilities. However, one could use efficient iteration methods [START_REF] Cavender | Quasi-stationary distributions of birth-and-death processes[END_REF][START_REF] Nåsell | Extinction and quasi-stationarity in the Verhulst logistic model[END_REF], approximations and asymptotic results. [START_REF] Kryscio | On the extinction of the S-I-S stochastic logistic epidemic[END_REF] summarized the first work on two popular approximations. The first one, denoted by p (0) , approximates the original birth and death process by assuming that μ 1 = 0; that is, the death rate of a population of size 1 is equal to 0 which amounts the original process with the origin 0 removed. In the context of the SIS epidemic model a second approximation, called p (1) , is obtained by assuming the existence of one permanently infected individual. The two approximations are birth and death processes without absorbing states, so their stationary distributions can be used to approximate the quasi-stationary distribution. In epidemiological models, the goodness of the approximation depends on the so-called transmission factor, which is denoted by R 0 . The transmission factor is interpreted in biological terms as the average number of contacts produced by one individual during the period of infectivity, when the individual is introduced into a completely susceptible population [START_REF] Allen | An Introduction to Stochastic Processes with Applications to Biology[END_REF]. The literature for the SIS model is particularly rich (see Subsection 3.1.1), but there also exists an important number of papers dealing with approximations and asymptotic results for other stochastic epidemic models (e.g. the Verhulst logistic model [START_REF] Nåsell | Extinction and quasi-stationarity in the Verhulst logistic model[END_REF], the SIR model with demography (Nåsell 1999b[START_REF] Nåsell | Extinction and quasi-stationarity in the Verhulst logistic model[END_REF]van Herwaarden and Grasman, 1995, etc.)).

We now concentrate on the case where S T is finite but reducible. This case has been investigated in a recent paper by [START_REF] Van Doorn | Survival in a quasi-death process[END_REF]. Some preliminary notation is needed. Suppose that S T consists of L communicating classes S k , for 1 ≤ k ≤ L. A partial order on {S k ; 1 ≤ k ≤ L} is defined by writing S i ≺ S j when class S i is accessible from S j . Let -α k be the (negative) eigenvalue with maximal real part of the sub-generator Q k corresponding to the states in S k . Then, the eigenvalue of Q S T with maximal real part is obtained as -α, where α = min 1≤k≤L α k . We also define I(α) = {k : α k = α} and a(α) = min I(α).

We can now summarize the main results in van Doorn and [START_REF] Van Doorn | Survival in a quasi-death process[END_REF]. If -α has a geometric multiplicity one, then the Markov chain has a unique quasi-stationary distribution {u i ; i ∈ S T } from which S a(α) is accessible. The jth component of {u i ; i ∈ S T } is positive (i.e., j ∈ S T + ) if and only if state j is accessible from S a(α) . A simple necessary and sufficient condition for establishing that -α has geometric multiplicity one is that {S k ; k ∈ I(α)} is linearly ordered, that is, S i ≺ S j ⇐⇒ i ≤ j, for all i, j ∈ I(α).

Following an example given by van Doorn and [START_REF] Van Doorn | Survival in a quasi-death process[END_REF], we notice that the quasi-stationary distribution of the pure death process on S = {0, 1, 2} with absorbing state 0 and μ i > 0, for i ∈ S T = {1, 2}, is given by

(u 1 , u 2 ) = (1, 0), if μ 1 ≤ μ 2 , μ 2 μ 1 , 1 -μ 2 μ 1 , if μ 2 < μ 1 . (2.3)
Suppose that X(0) = 2 and μ 1 = μ 2 = μ, then formula (2.3) says that the process will be in state i = 1 with probability 1, if after a long time absorption has not yet occurred. In other words, the set of states having positive mass is reduced to S T + = {1}. However, the absorption certainly occurs as far as t → ∞, and the expected amount of time that the processes spends in state i ∈ S T = {1, 2} is 1/μ, given that X(0) = 2. Thus, one could appeal to a natural frequency interpretation of the concept of probability to conclude that the even distribution that gives mass 0.5 to each state i ∈ S T = {1, 2} provides the natural measure of the process behavior at an arbitrary time before absorption. However, if one wishes to know the state distribution after a very long time, given that the process is known to be in S T , then the quasi-stationary distribution is the natural solution.

Comparing the dynamics of the pure death model versus most epidemic models, one may conclude that they are so different. It is certainly true, but our interest in the pure death model is just that it provides an extremely simple scenario where the quasi-stationary distribution and the RE-distribution exhibit drastically different solutions. The first, but not the unique, key to understand this peculiar fact is that S T is reducible.

The concept of quasi-stationarity becomes more complicated when S T is infinite. Then, a natural question of interest is to find finite approximations. Let us assume that S T is irreducible. Then, under appropriate conditions [START_REF] Breyer | Approximations of quasi-stationary distributions for Markov chains[END_REF] there exists an increasing sequence of finite irreducible truncations whose quasi-stationary distributions are used to approximate the corresponding one of the original infinite process. The parallel discrete-time results can be found in [START_REF] Seneta | Non-Negative Matrices and Markov Chains[END_REF].

We now turn our attention to the reformulation of the RE-distribution. Let T j be the time that process {X(t); t ≥ 0} spends in state j ∈ S T before absorption. Then, given that X(0) = i ∈ S T , we define

P i (j) = E i [T j ] E i [T ]
, i, j ∈ S T .

(2.4)

In the light of (2.4), we notice that the RE-distribution always exists whether S T is finite or infinite, provided that the expected time to absorption

E i [T ] = E [T | X(0) = i ] is finite.
Moreover, it should be pointed out that the REdistribution assigns positive probability to all state j accessible from the initial state i. Compared to the quasi-stationary distribution, this is an important difference.

The RE-distribution has a meaningful probabilistic interpretation. Let us construct the ideal replicated model obtained by assuming that at each extinction the biological model restarts in the same initial state i ∈ S T . This replicated model is a regenerative process so its limiting distribution equals the expected amount of time the process spends in the state j during one regeneration cycle divided by the expected length of one cycle, see Theorem 8.26 in [START_REF] Kulkarni | Modeling and Analysis of Stochastic Systems[END_REF]. In Subsection 3.1.1, we will explore in more detail the relationship between the quasi-stationary and ratio of expectation distributions, and the model that immediately restarts in state i.

The ratio of means distribution studied in the early literature is now obtained as the following unconditional version of the RE-distribution (2.4):

P (j) = i∈S T a i E i [T j ] i∈ST a i E i [T ]
, j ∈ S T .

(2.5)

In the next section we deal with a variant of P (j) in (2.5) obtained by replacing the initial distribution {a i ; i ∈ S T } by weights {ω i ; i ∈ S T }, which do not preserve the interpretation as initial probabilities at time t = 0. Darroch andSeneta (1965, 1967) do not pursue the use of {P (j); j ∈ S T } further because it depends on the initial distribution {a i ; i ∈ S T }. More concretely, they clearly stated that the reason for the dependence of the ratio of means on the initial distribution is that the absorption time is not long enough for the dependence on the initial distribution to wear off. Darroch and Seneta (1965, Section 4) also explain that the quasi-stationary distributions do not depend on the initial distribution and are derived, roughly speaking, by considering only those realizations in which the absorption time is long. Since then, the ratio of means distribution has been largely ignored.

There is no doubt that the dependence on the initial distribution is a significant conceptual difference between both approaches. In population biology, it is often known that a certain population has been evolving for a long time, and that it has not reached the extinction yet. Furthermore, it may be very difficult, or even impossible, to know the exact initial distribution. From this perspective, the quasi-stationary distribution gives a very satisfactory approximation of the system state before absorption (i.e., the first objective mentioned in the abstract). However, we think that, whenever it can be possible, the effect of the initial distribution should be taken into account. We note that formula (2.4) allows the use of the RE-distribution not only when the time to absorption is large (i.e., the region R 0 > 1 in epidemic models) but also in cases where absorption is reached soon. If the absorption time is not long enough, then the initial distribution influences the current system state. In this case, the efforts to deal with the RE-distribution are welcome. For example, an outbreak of head lice (pediculosis capitis) in a school should end within a not so long time horizon, and the school administrators should have some knowledge about the number of pupils infected when the outbreak was detected. Despite of whether it can be observed or not, it seems reasonable to assume that an outbreak of many epidemics starts with a single infective who introduces the infection into the population. Bayesian methods [START_REF] Clancy | Bayesian estimation of the basic reproduction number in stochastic epidemic models[END_REF] and the maximum entropy principle can be helpful for estimating initial distributions. On the other hand, in those situations where the system parameters lead to long absorption times, the effect of the initial state X(0) = i is weak. The numerical examples in Section 3 illustrate that the quasi-stationary distribution and the RE-distribution are in fact close when the absorption time is long. In this important case, the RE-distribution can be used as a simple approximation to the quasi-stationary distribution (i.e., the second objective in this paper).

To conclude this section, we next discuss how to compute the RE-distribution.

Obviously, E i [T ] = ∞ if P {T < ∞ | X(0) = i } < 1. In the case P {T < ∞ | X(0) = i } = 1 (i.e.
, absorption occurs with probability 1), the expectations {E i [T ] ; i ∈ S T } are given by the smallest non-negative solution to the system of linear equations (see Theorem 6.19 in [START_REF] Kulkarni | Modeling and Analysis of Stochastic Systems[END_REF])

E i [T ] = 1 q i + j∈S T j =i q ij q i E j [T ] , i ∈ S T . (2.6)
Then, by expressing (2.6) in matrix form, we obtain

Q S T m = -e,
where m and e are column vectors with dimension the cardinality of the subset S T . The entries of m are E i [T ] , whereas all entries of e are all equal to 1. For a fixed j ∈ S T , we employ a first-step argument to find that

E i [T j ] = δ ij q i + k∈S T k =i q ik q i E k [T j ] , i ∈ S T , (2.7)
where δ ij denotes Kronecker's delta which is defined as follows:

δ ij = 1, if i = j, 0, otherwise.
Thus, the matrix form version of equation (2.7) is

Q S T m j = -e j , (2.8)
where m j is now the column vector with entries E i [T j ] , for i ∈ S T , and e j is a column vector such that all entries are equal to 0, except the jth which is equal to 1.

Applications to stochastic population models

In Section 2 we introduced the RE-distribution in the context of an absorbing continuous-time Markov chain with countable state space. In this section we provide applications to some selected population models: the classical SIS epidemic model, a linear growth model, the classical SIR model and a finite SIR model with demography. A comparison with the quasi-stationary distribution is done. Our results include a few theoretical findings and numerical experiments. The main results for the classical SIS and SIR epidemic models can be found in several textbooks [START_REF] Andersson | Stochastic Epidemic Models and Their Statistical Analysis[END_REF][START_REF] Allen | An Introduction to Stochastic Processes with Applications to Biology[END_REF][START_REF] Bailey | The Elements of Stochastic Processes with Applications to the Natural Sciences[END_REF][START_REF] Daley | Epidemic Modelling: An Introduction[END_REF]. For discussions about variants and generalizations, please consult some recent papers [START_REF] Allen | An Introduction to Stochastic Epidemic Models[END_REF][START_REF] Clancy | A stochastic SIS infection model incorporating indirect transmission[END_REF][START_REF] Fouchet | Impact of myxomatosis in relation to local persistence in wild rabbit populations: The role of waning immunity and the reproductive period[END_REF][START_REF] Lindholm | On the time to extinction for a two-type version of Bartlett's epidemic model[END_REF][START_REF] Xu | Stochastic model of an influenza epidemic with drug resistance[END_REF] and the references therein.

General birth and death processes

We consider a birth and death process {X(t); t ≥ 0} with state space S = {0, 1, ..., N }. Both the finite case (N < ∞) and the case of infinite state space (N = ∞) can be studied. The infinitesimal transition rate from state i to state i + 1 is denoted by λ i > 0 (birth rate), for 1 ≤ i ≤ N -1, whereas the transition rate from state i to state i -1, for 1 ≤ i ≤ N , is denoted by μ i > 0 (death rate). We assume that the origin is an absorbing state, so λ 0 = 0. Moreover, if N < ∞ then we take λ N = 0 to be consistent with the assumption that S is finite. Figure 1 shows the state space and transitions.

Fig. 1. States and transitions of the birth and death model

• • • • • 0 1 2 N -1 N μ 1 - λ1 μ 2 - λ2 • • • μ N -1 - λN-1 μ N For each fixed j ∈ {1, ..., N }, equation (2.7) yields E 0 [T j ] = 0, (3.1) (λ i + μ i )E i [T j ] = μ i E i-1 [T j ] + λ i E i+1 [T j ] , i = j, 1 ≤ i ≤ N, (3.2) (λ j + μ j )E j [T j ] = μ j E j-1 [T j ] + λ j E j+1 [T j ] + 1. (3.3)
Using methods of finite difference equations, we find that the solution of the system (3.1)-(3.3) is given by

E i [T j ] = 1 μ j min(i,j) k=1 j-1 n=k λ n μ n , 1 ≤ i, j ≤ N. (3.4) Since (T | X(0) = i ) = N j=1 (T j | X(0) = i ), summing (3.4
) over j we obtain

P i (j) = E i [T j ] E i [T ] = 1 μ j min(i,j) k=1 j-1 n=k λn μ n N j=1 1 μ j min(i,j) k=1 j-1 n=k λn μ n , 1 ≤ i, j ≤ N. (3.5)
We emphasize that explicit formula (3.5) for the RE-distribution applies both for N < ∞ and N = ∞. Furthermore, after simple algebraic manipulations, it is easy to observe that the expression given in the denominator for E i [T ] is in agreement with the well-known expression for the expected time until extinction, see formula (5.7) in [START_REF] Norden | On the distribution of the time to extinction in the stochastic logistic population model[END_REF] and formula (6.22) in [START_REF] Allen | An Introduction to Stochastic Processes with Applications to Biology[END_REF].

We conclude this subsection by noticing that the whole probability distribution of the variables T j | X(0) = i can be easily determined in terms of a system of linear equations for the Laplace transforms

E [exp {-sT j } | X(0) = i ] .
Once the Laplace transforms have been computed, the corresponding density functions, or the survival probabilities P i {T j > x}, can be obtained by numerical inversion methods [START_REF] Cohen | Numerical Methods for Laplace Transform Inversion[END_REF]. Equations for higher moments E i T k j follow by direct differentiation of the Laplace equations.

The classical SIS model

The general formulation given in Subsection 3.1 includes as particular cases many stochastic population models. In particular, for the case N < ∞, we mention the stochastic version of the Verhulst logistic model and the stochastic SIS epidemic model [START_REF] Allen | An Introduction to Stochastic Processes with Applications to Biology[END_REF]. The latter model corresponds to the following specification:

λ i = β N i(N -i), 0 ≤ i ≤ N, μ i = (γ + b)i, 1 ≤ i ≤ N,
where β, γ and b respectively denote the contact, recovery and birth-death rates.

We assume that b = 0. Then, the transmission factor is R 0 = β/γ. The quasi-stationary distribution of the SIS model cannot be obtained explicitly. However, the probabilities {u i ; 1 ≤ i ≤ N } can be computed by using the recursive procedure proposed by [START_REF] Cavender | Quasi-stationary distributions of birth-and-death processes[END_REF]. Two methods for implementing the recursion are summarized by [START_REF] Nåsell | Extinction and quasi-stationarity in the Verhulst logistic model[END_REF]. In addition, the two approximations p (0) and p (1) described in Section 2 are often useful. The quality of the approximations is discussed with respect to the three parameter regions where R 0 > 1 (R 0 is distinctly larger than 1), R 0 < 1 (R 0 is distinctly smaller than 1) and R 0 is close to 1 (transient region), as N → ∞. Following Nåsell (1999a), we notice that in the regions R 0 > 1 and R 0 < 1 the quasistationary distribution is well approximated by p (0) and p (1) , respectively. Nåsell (1999a) derives a more sophisticated approximation for the intricate transient case. From a practical point of view, the transmission factor R 0 is parametrized as R 0 = 1 + ρ/ √ N. Then, Nåsell proposes to use the somewhat arbitrary value ρ = 3 to determine the boundaries of the three regions.

On the other hand, in the case of the SIS model, the RE-distribution (3.5) reduces to

P i (j) = 1 j min(i,j) k=1 R0 N j-k (N -k)! (N -j)! N j=1 1 j min(i,j) k=1 R0 N j-k (N -k)! (N -j)! , 1 ≤ i, j ≤ N. (3.6)
At a first glance, we easily derive from (3.6) the epidemic behavior as far as the transmission factor increases but the population size remains constant. Then, we expect long absorption times while the mass probability is concentrated on the state N . This yields lim

R0→∞ E i [T ] = ∞, 1 ≤ i ≤ N, lim R0→∞ P i (j) = δ jN , 1 ≤ i, j ≤ N, lim R0→∞ N j=1 jP i (j) = N, 1 ≤ i ≤ N.
We next summarize a set of stochastic ordering results which are useful to compare the quasi-stationary distribution, the RE-distribution and the birth and death approximations p (0) and p (1) . For notational convenience, in the sequel we employ the bold notation u, P i , for 1 ≤ i ≤ N, p (0) and p (1) to denote the row vectors containing the probabilities of the corresponding distribution. The symbol ≤ st will denote the usual stochastic order relation with respect to the distribution function.

First, we notice that

P 1 = p (0) , (3.7)
as it was to be expected from the probabilistic interpretation of both distributions.

We now observe that

p (0) ≤ st u ≤ st p (1) , (3.8) P 1 ≤ st u ≤ st P N ,
(3.9)

P i ≤ st P i , 1 ≤ i ≤ i ≤ N.
(3.10)

Relation p (0) ≤ st u was proved by [START_REF] Cavender | Quasi-stationary distributions of birth-and-death processes[END_REF] for the general birth and death process. For the SIS model, [START_REF] Kryscio | On the extinction of the S-I-S stochastic logistic epidemic[END_REF] observed empirically that p (1) is an upper bound of u. This conjecture has been proved by [START_REF] Clancy | A note on quasi-stationary distributions of birth-death processes and the SIS logistic epidemic[END_REF] with the help of an auxiliary process that, whenever the system reaches the absorption, it automatically restarts in state j, for 1 ≤ j ≤ N , with probability ν j . In particular, by assigning unit mass to a fixed state i (i.e., we take ν j = δ ij ), we obtain formula (3.9); see Corollary 1 in [START_REF] Clancy | A note on quasi-stationary distributions of birth-death processes and the SIS logistic epidemic[END_REF]. The result in (3.10) follows from Theorem 1 in that paper. It turns out that formulas (3.9) and (3.10) hold for the general birth and death process.

At this point, we also observe that the RE-distribution P i is in fact equal to the stationary distribution of the auxiliary process with ν j = δ ij , for 1 ≤ j ≤ N . In this way, we complete the circle arriving again at the probabilistic interpretation for P i already given in Section 2, where we appealed to the relationship with the regenerative processes.

Formulas (3.8) and (3.9) show that p (1) and P N are upper bounds of the quasi-stationary distribution u. Thus, it is appropriate to study the relation between them. To this end, for N fixed, we find that lim

R0→∞ l j=1 p (1) j l j=1 P N (j) = l N ≤ 1, 1 ≤ l ≤ N, lim R0→0 l j=1 p (1) j l j=1 P N (j) = N j=1 1 j l j=1 1 j ≥ 1, 1 ≤ l ≤ N,
which amounts to the following relations: (1) , for N fixed and R 0 sufficiently large, (3.11) p (1) ≤ st P N , for N fixed and R 0 sufficiently small. (3.12)

P N ≤ st p
From (3.9) we see that the extreme RE-distributions P 1 and P N are respectively lower and upper bounds of the quasi-stationary distribution u. Moreover, formula (3.10) shows that the RE-distributions are stochastically monotone with respect to the initial state. Therefore, one may conjecture that either any appropriate choice of i or the consideration of a mixture P = N i=1 ω i P i would lead to some improvement, if one wishes to consider the RE-distribution as an approximation of the quasi-stationary distribution. In order to explore numerically this possibility, we next suggest two RE-approximations.

A first approximation, denoted by P, approximates u by a single appropriate RE-distribution P i . The choice of the appropriate index i depends on the three parameter regions described earlier. When R 0 is distinctly larger than 1, we use the good approximation p (0) to determine the starting state as i = E[p (0) ] , where x is the nearest integer function defined as the closest integer to x. In the parameter region where R 0 is distinctly smaller than 1, we use p (1) to get i = E[p (1) ] . In the more involved transient region, we start from a rough but simple interpolation approximation p = R0 1+R0 p (0) + 1 1+R0 p (1) . Large values of the first weight give support to long extinction times, in agreement with the use of the approximation p (0) . Then, we take i

= E[p] .
As a second approximation, we propose to employ a mixture of the REdistributions, where the weights ω i are given by the first RE-approximation P. In this way, we construct the distribution P, which assigns probabilities Some numerical results are summarized in Tables 1 and2. All codes have been done using Fortran. In Table 1 we compare six distributions, namely the quasi-stationary distribution u, the birth and death approximation p (1) , the extreme RE-distributions P 1 and P N , and the two RE-approximations P and P. From (3.7), we recall that p (0) = P 1 . We choose the population size as N = 100 and normalize the recovery rate to be equal to one. Then, we take R 0 = 0.5, 0.9, 1.0, 1.3, 1, 5 and 2.0. The choice here covers the three parameter regions. More concretely, we interpret that the cases R 0 = 0.9 and 1.0 belong to the transient region, while R 0 = 1.3 is somewhat in the boundary with the region where R 0 is distinctly larger than 1. Each entry in the table gives, from top to bottom, the mean value, the standard deviation, the first probability and the mode, for the corresponding distribution. For example, in the case of u the notation of these characteristics is E[u], σ(u), u 1 and u m . Moreover, the last row gives the expected time to extinction from the states 1 and N. The numbers in bold indicate what is the closest distribution to the quasi-stationary one, with respect to the corresponding characteristic.

P j = N i=1 P i P i (j), for 1 ≤ j ≤ N. R 0 0.5 0.9 1.0 1.3 1.5 2.0 E[u] σ(u) u 1 u m 1.
As a supplement to Table 1, we notice that R u decreases from the value 0.95899, for R 0 = 0.5, to the value 1.79 × 10 -8 , corresponding to R 0 = 2.0. In fact, R u = 0.99999, for R 0 = 10 -9 . Thus, we conclude that the quasi-stationary regime in the SIS model is always reached faster than the eventual absorption.

A summary of the observations inferred from Table 1 is given below.

• It is clear that p (1) gives the best approximation in the region where R 0 is distinctly smaller than 1.

• Regarding the transient region, we remark the interest of an asymptotic approximation due to Nåsell (1999a), which is clearly better than p (0) and p (1) . For the sake of simplicity, we have avoided the implementation of that approximation. Among the distributions under comparison here, we notice that the RE-approximations P and P are satisfactorily close to u in the cases R 0 = 1.0 and 1.3. Thus, it seems reasonable to consider the use of P and P in the interval 0 ≤ ρ ≤ 3, specially when N is moderate or small and there exist doubt about the applicability of asymptotic results.

• It is remarkable that all the distributions, except p (1) , are satisfactory approximations of u in the important practical region where R 0 is distinctly larger than 1. However, the best comparative indicators correspond to the distribution P. The differences among the distributions are more notorious as far as ρ tends to 3. In this case, it would be possible to use either P or P.

• Combining the stochastic ordering relationships (3.8)-(3.12) and the known asymptotic results [START_REF] Nåsell | The quasi-stationary distribution of the closed endemic SIS model[END_REF](Nåsell , 1999a)), it is easy to conclude that the mean values of the six distributions under study are asymptotically equivalent to N (R 0 -1)/R 0 , in the region R 0 > 1. The numerical results in the table support this observation.

• Note that formula (7.3) in [START_REF] Nåsell | The quasi-stationary distribution of the closed endemic SIS model[END_REF] gives the approximation

u 1 = √ N (R 0 -1) 2 ϕ(β 1 )/R 0 , where ϕ(β 1 ) = 1 √ 2π e -N log R0- R 0 -1 R 0
. When R 0 = 1.3, we have u 1 = 0.0169, while u 1 = 8.1583 × 10 -9 for R 0 = 2.0. An explanation of the superiority of the RE-approximations over u 1 is the following. From the approximation (8.4) for E N [T ] in [START_REF] Nåsell | The quasi-stationary distribution of the closed endemic SIS model[END_REF], we find that P N (1) can be approximated by

P N (1) = √ NR N -1 0 (R 0 -1) 2 ϕ(β 1 )/(R N 0 -1)
. Now we observe that u 1 < P N (1) but this disagrees with formula (3.9), which implies the exact relation P N (1) ≤ u 1 . The interest in the first probability u 1 comes from its relation with the expected time to extinction when the initial distribution is u. Denoting the latter as E[T u ], we recall that both quantities are related through the formula E[T u ] = (γu 1 ) -1 .

In Table 2, we employ the maximum pointwise distance to compare the distributions under consideration versus the quasi-stationary distribution. For example, in the case of p (1) , the distance is defined by p (1)u = max 1≤j≤N p (1) ju j . The entry marked in bold is again associated with the distribution providing the minimum distance and, consequently, the best global approximation. R 0 0.5 0.9 1.0 1. Table 2. Distributions distances with respect to u

In the light of the table, we see that the best approximation starts being p (1) , in the region R 0 < 1, but passes by P, in the transient region, and finally becomes P, in the region R 0 > 1. It should be pointed out the significant improvement obtained in the region R 0 > 1 when we use P.

A linear growth model

As an example of the infinite case N = ∞, we next consider the linear growth model with transition rates

λ i = λi, i ≥ 0, μ i = μi, i ≥ 1.
Linear transition rates arise often in combination with catastrophes, see for instance [START_REF] Lee | The density of the extinction probability of a time homogeneous linear birth and death process under the influence of randomly occurring disasters[END_REF]. Moreover, the linear growth model is one of the few cases where both the quasi-stationary distribution and the RE-distribution can be determined explicitly [START_REF] Nåsell | Extinction and quasi-stationarity in the Verhulst logistic model[END_REF]. This fact facilitates a rapid comparison between both distributions.

We assume that R 0 = λ/μ < 1. Then, the quasi-stationary distribution follows the geometric law:

u j = (1 -R 0 )R j-1 0 , j ≥ 1.
On the other hand, formula (3.5) for the RE-distribution yields

P i (j) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ (1-R j 0 )j -1 (1-δ1i) i-1 j=1 1-R j-i 0 j - 1-R i 0 R i 0 ln(1-R0) , if 1 ≤ j < i, R j-i 0 (1-R i 0 )j -1 (1-δ1i) i-1 j=1 1-R j-i 0 j - 1-R i 0 R i 0 ln(1-R0) , if 1 ≤ i ≤ j.
It is easy to show that both distributions have a unique mode at the point u m = P m i = 1, for i ≥ 1, whereas the expected values are given by

∞ j=1 ju j = 1 1 -R 0 , ∞ j=1 jP i (j) = i ⎛ ⎝ (1 -δ 1i ) i-1 j=1 1 -R j-i 0 j - 1 -R i 0 R i 0 ln(1 -R 0 ) ⎞ ⎠ -1 , i ≥ 1.
We are interested in the limiting behavior as R 0 → 1. Then, we obtain the following limiting results:

lim R0→1 E i [T ] = ∞, i ≥ 1, (3.13) lim R0→1 ∞ j=1 ju j = ∞, (3.14) lim R0→1 ∞ j=1 jP i (j) = ∞, i ≥ 1, (3.15) lim R0→1 ∞ j=1 jP i (j) ∞ j=1 jP i (j) = 1, 1 ≤ i, i ≤ N, (3.16) lim R0→1 ∞ j=1 jP i (j) ∞ j=1 ju j = 0, i ≥ 1.
(3.17) Formula (3.13) says that the expected absorption time becomes arbitrarily long, as far as R 0 → 1, so the comparison is meaningful. In agreement with the intuition, formulas (3.14) and (3.15) show that the mean value of the two distributions tends to infinity. From (3.15), we also notice that the influence of the initial state i vanishes when the time to absorption is long. Formula (3.16) says that the mean values of the RE-distribution obtained for two different initial states are asymptotically equivalent, when R 0 → 1. Finally, the limiting ratio in (3.17) illustrates the important fact that the mean value of the quasistationary distribution converges to infinity faster than the mean value of the RE-distribution. In this average sense, a similar limiting behavior of both distributions cannot be concluded.

SIR epidemic models

The SIR formalism categorizes individuals within a closed population as susceptible, infective and removed individuals. The main difference between SIS and SIR models is that, in the SIR model, individuals recover and develop permanent immunity. We are concerned with the classical SIR model and with a finite SIR model with demography. The classical SIR model refers to the standard formulation as it can be found in many textbooks [START_REF] Andersson | Stochastic Epidemic Models and Their Statistical Analysis[END_REF][START_REF] Allen | An Introduction to Stochastic Processes with Applications to Biology[END_REF][START_REF] Bailey | The Elements of Stochastic Processes with Applications to the Natural Sciences[END_REF][START_REF] Daley | Epidemic Modelling: An Introduction[END_REF]. The introduction of demographic forces in the epidemiological formalism was already proposed by [START_REF] Bartlett | Deterministic and stochastic models for recurrent epidemics[END_REF]. Slight variants of Bartlett's model lead to what is known as the SIR model with demography [START_REF] Van Herwaarden | Stochastic epidemics: major outbreaks and the duration of the endemic period[END_REF]Nåsell, 1999b). These models and their generalizations (see e.g. [START_REF] Nåsell | Stochastic models of some endemic infections[END_REF] are based on Markov chains with infinite state space. In Subsection 3.2.2, we introduce a finite SIR model with demography. The finitess of the state space makes the model enough tractable, which is needed to deal with our objective of comparing the quasi-stationary distribution and the ratio of expectations distribution.

The existence of a finite state space, which can be partitioned into two classes of absorbing and transient states, is a common feature for the two SIR models under study, and even for the classical SIS model. However, the quasistationary behavior of these models is completely different due to the different dynamics of the transitions among the transient states.

The classical SIR model

Let {(X(t), Y (t)); t ≥ 0} be the bidimensional continuous-time Markov chain describing the SIR epidemic model. At time t, the population consists of X(t) infectives, Y (t) susceptibles and Z(t) = N -X(t) -Y (t) immune individuals, where N is the constant population size. The initial condition is (X(0), Y (0)) = (m, n). When in state (i, j), for i ≥ 1, the population state moves either to (i + 1, j -1) at rate λ ij (λ i0 = 0) due to an infection, or to (i -1, j) at rate μ i (μ 0 = 0) due to the removal of an infective. The state space of the SIR epidemic model is S = {(i, j); 0 ≤ i ≤ m + n, 0 ≤ j ≤ min{n, m + n -i}}. We notice that states (0, j), for 0 ≤ j ≤ n, are absorbing states, so it is reasonable to assume that m ≥ 1. Fig. 2 illustrates the transitions among the states for the case (m, n) = (3, 3).

0 1 2 3 Y(t) 0 1 2 3 4 5 6 X(t)

Fig. 2. States and transitions of the SIR epidemic model

• • • • • • • • • • • • • • • • • • • • • • μ 1 μ 2 μ 3 μ 4 μ 5 μ 6 μ 1 μ 2 μ 3 μ 4 μ 5 μ 1 μ 2 μ 3 μ 4 μ 1 μ 2 μ 3 @ @ R @ @
λ11 @ @ R @ @ λ21 @ @ R @ @ λ31 @ @ R @ @ λ41 @ @ R @ @ λ51 @ @ R @ @ λ12 @ @ R @ @ λ22 @ @ R @ @ λ32 @ @ R @ @ λ42 @ @ R @ @ λ13 @ @ R @ @ λ23 @ @ R @ @ λ33 A typical choice for the transition rates is [START_REF] Allen | An Introduction to Stochastic Processes with Applications to Biology[END_REF])

λ ij = β N ij, (i, j) ∈ S, μ i = γi, 0 ≤ i ≤ m + n,
where β and γ denote the contact and the recovery rates.

From the results summarized in Section 2 it follows easily that the quasistationary distribution only assigns positive probabilities to the states (i, 0) with 1 ≤ i ≤ i * , where i * is the minimum i, for 1 ≤ i ≤ m + n, such that μ i = min 1≤k≤m+n μ k . In particular, if μ i = γi, for 1 ≤ i ≤ m + n, then i * = 1 and the quasi-stationary distribution has the trivial form

u (i,j) = δ (1,0)(i,j) , 0 ≤ j ≤ n, 1 ≤ i ≤ m + n -j,
that is, all the probability mass is concentrated at the state (1, 0).

On the other hand, the RE-distribution of the SIR epidemic model has the form

P (m,n) (i, j) = E (m,n) T (i,j) E (m,n) [T ] , 0 ≤ j ≤ n, 1 ≤ i ≤ m + n -j. (3.18)
Define A ij to be the probability of reaching the state (i, j) ∈ S starting from (m, n) before the extinction occurs. Considering that each state (i, j) can be visited at most one time, we write formula (3.18) as follows:

P (m,n) (i, j) = Aij λij +μ i n j=0 m+n-j i=1 Aij λij +μ i , 0 ≤ j ≤ n, 1 ≤ i ≤ m + n -j.
Next we summarize the recursive scheme for the computation of the probabilities A ij [START_REF] Neuts | An algorithmic study of S-I-R stochastic epidemic models[END_REF].

Step 1. Set A ij = δ (i,j)(m,n) .

Step 2. For i = m -1, m -2, ..., 0 calculate

A in = A i+1,n μ i+1 λ i+1,n + μ i+1 .
Step 3. For k = 1, calculate 3.a.

A m+k,n-k = A m+k-1,n-k+1 λ m+k-1,n-k+1 λ m+k-1,n-k+1 + μ m+k-1 . 3.b. For i = m + k -1, m + k -2, ..., 2 compute A i,n-k = A i+1,n-k μ i+1 λ i+1,n-k + μ i+1 + A i-1,n-k+1 λ i-1,n-k+1 λ i-1,n-k+1 + μ i-1 .
3.c For i = 0, 1 calculate

A i,n-k = A i+1,n-k μ i+1 λ i+1,n-k + μ i+1 .
Step 4. Set k = k + 1. If k ≤ n go to Step 3.a.

We notice that the above recursive scheme also provides a method for computing the final size of the epidemic. In fact, A 0j , for 0 ≤ j ≤ n, gives the probability that there are j susceptible individuals at the extinction time. Thus, the final size of the epidemic is Nj.

At this point, we would like to note that the RE-distribution gives positive probability to all transient states S T = S -{(0, j); 0 ≤ j ≤ n}. In contrast, the typical choice of the SIR gives all the quasi-stationary probability mass to (1, 0). As in the example for the simple death process showed in Section 2, the explanation is related to the fact that S T is reducible, but also to the fact that state (1, 0) has the lowest exit rate (i.e., min (i,j)∈S T (λ ij + μ i ) = λ 10 + μ 1 = γ). Moreover, once a state has been visited, the system leaves it forever. In this sense, the SIR formalism is more appropriate to model acute infections including chickenpox, measles, mumps and rubella [START_REF] Keeling | Modeling Infectious Diseases in Humans and Animals[END_REF]), where in a relatively rapid period the extinction time is reached. The ratio R u gives support to the same idea. It is easy to prove that the eigenvalues of Q S T in the SIR model are given by -(λ ij + μ i ), for 0 ≤ j ≤ n and 1 ≤ i ≤ m + nj. Then, we readily find that R u ≥ 1. More concretely, we have

R u = 1, if γ ≤ β/N, 2γ β/N+γ , if γ > β/N.
In the light of the above comments, it should be erroneous to conclude that the concentration of the mass in a single point is a drawback of the quasistationarity. The correct interpretation is that it is not meaningful to study the quasi-stationary distribution in this case, because probably the time to absorption is not long enough. In contrast, the RE-distribution can be used despite how long the absorption time is, so it provides a natural approach to measure the behavior of the SIR model before absorption. Fortunately, we have more than one option to cover our basic objective.

A finite SIR model with demography

In the SIR model with demography (Nåsell, 1999b), the population is subject to immigration and deaths. It is usual to assume a constant birth rate μN and a linear death rate μj, proportional to the number j of susceptibles. Note that N is not the total population size. In fact, the population size varies along the time, but N can be viewed as the mean equilibrium population size when the infection has been eradicated. The state space of the Markov chain modelling the SIR model with demography is S = N × N and the states {(0, j); j ≥ 0} form the absorbing set S A . In this context, the quasi-stationary distribution cannot be determined neither explicit nor recursively. The study of the time to extinction is also intricate, even to prove almost sure extinction is an involved problem. The existing work (Nåsell, 1999b[START_REF] Nåsell | Stochastic models of some endemic infections[END_REF] is mainly concerned with the determination of asymptotic approximations of the quasi-stationary distribution and the time to extinction, when the ratio R 0 is distinctly larger than 1. We next consider alternative demographic rates leading to a Markov chain with finite state space S = {(i, j); 0 ≤ i+j ≤ N }. With this state space, the comparative study of the quasi-stationary and the ratio of expectations distributions can be numerically implemented.

We consider a bidimensional process {(X(t), Y (t)); t ≥ 0}, where X(t) denotes the number of infected individuals, while Y (t) represents the number of susceptibles at time t. When in state (i, j), the dynamics of the system is defined by the following four transitions. Due to an infection the population size moves to (i + 1, j -1) at rate λ ij = βij/N , while the removal of an infective causes a transition to (i -1, j) at rate μ i = (γ + μ + θ)i. Here β and γ represent the contact and the recovery rate, as in the classical SIR model. The rates μ and θ correspond to natural death and death due to the infection, respectively. Due to the demographic forces the system moves either to state (i, j + 1) at rate

α ij = ( λ N (i + j) + ξ)(N -i -j)
, when a birth takes places, or to (i, j -1) at rate δ j = μj, due to the death of a susceptible. The rate α ij can be decomposed into two contributions. The first term λ N (i + j)(Nij) shows that the internal births follow a logistic growth. With this rate, the population cannot growth beyond the carrying capacity N . On the other hand, the term ξ(Nij) gives the rate of the external (immigration) births. We assume that the arrival of immigrants decreases as far as the population tends to its ideal carrying capacity. The transition among states are represented in Figure 3, for the simple case N = 3.

• • • • • • • • • • • • • • • μ 1 μ 2 μ 3 μ 4 μ 1 μ 2 μ 3 μ 1 μ 2 μ 1
@ @ @ R @ @ @ λ11 @ @ @ R @ @ @ λ21 @ @ @ R @ @ @ λ31 @ @ @ R @ @ @ λ12 @ @ @ R @ @ @ λ22 @ @ @ R @ @ @ λ13 6 α00 ? δ1

6 α01 ? δ2 6 α02 ? δ3 6 α03 ? δ4 6 α10 ? δ1 6 α11 ? δ2 6 α12 ? δ3 6 α20 ? δ1 6 α21 ? δ2 6 α30 ? δ1
As a partially related work, we mention the papers by [START_REF] Swift | A logistic birth-death-immigration-emigration process[END_REF] and [START_REF] Allen | Comparison of deterministic and stochastic SIS and SIR models in discrete time[END_REF]. The former deals with a model operating in the presence of logistic growth with immigration. In the latter, an SIR model with finite triangular state space is considered. In that paper, each death is accompanied by a birth so that the population size remains constant.

In Table 3, we illustrate the main characteristics of the quasi-stationary probabilities u ij , for (i, j) ∈ S T = {(i, j); 1 ≤ i+j ≤ N }, which can be computed from the fact that the quasi-stationary distribution u is the left eigenvector of Q S T with maximal real part. The finite SIR model with demography has seven system parameters: N, β, γ, μ, θ, λ and ξ. In our numerical experiment, we fix N = 100, γ = 499/500, μ = 1/500, θ = 0 and ξ = 0.1. The choice of γ and μ normalizes γ + μ as 1, and reflects that the expected life is much longer than the expected time to recovery. The internal birth rate λ is chosen to be multiple of ξ; that is, we consider λ = 0.1, 0.5, 1.0, 5.0 and 10.0. Finally, the transmission factor R 0 = β/(γ + μ) = β takes values R 0 = 0.5, 1.0, 1.25, 1.5 and 2.0. For each pair (R 0 , λ), we compute the marginal expected values of the number of infectives, E[u i. ], and the number of susceptibles, E[u .j ], the standard deviation of the number of infectives, σ(u i. ), and the mode u m ij (i.e., the bidimensional pair (i, j) with highest quasi-stationary mass). An examination of the table reveals that E[u i. ] is an increasing function of R 0 and λ. The mean value E[u .j ] decreases as a function of R 0 and it has a mode as a function of λ, in the region R 0 > 1. The standard deviation σ(u i. ) and the mode u m ij exhibit higher magnitudes in the southeastern corner of the table.

E[u

i. ] E[u .j ] σ(u i. ) u m ij R 0 = 0.5 R 0 = 1.0 R 0 = 1.25 R 0 = 1.5 R 0 = 2.0 λ = 0.
In Table 4 we are concerned with the P distribution obtained by starting the RE-distribution P (m,n) (i, j) at the point ( m, n), such that m = E[u i. ] and n = E[u .j ] . The P distribution can be computed by solving numerically the systems of linear equations that govern the expectations E (m,n) T (i,j) and E (m,n) [T ] . The block tridiagonal structure of the infinitesimal generator might be exploited to this end. From our numerical experiments, we now observe that R u < 1 for all choices of the pair (R 0 , λ) in the tables, except for the case (R 0 , λ) = (0.5, 0.1), where R u = 1.41443. As a result, the comparison between P and u is meaningful in a broad domain of the considered parameters.

By comparing the entries in Tables 3 and4, we may conclude that P provides an acceptable approximation of u in the whole range of the pair (R 0 , λ). The quality of the approximation is improved when R 0 and λ simultaneously increase and, consequently, the time to extinction becomes longer. In fact, both distributions present identical values in the case (R 0 , λ) = (2.0, 10.0). Numerical experiments, not reported here, show that the characteristics of the P distribution are almost insensitive with respect to the initial state, as far as R 0 and λ increase.

Concluding comments

Our aim in this paper is to provide a comparative analysis of the quasistationary and the ratio of expectation distributions. The motivation for this study comes from the need of understanding the behavior of a biological stochastic system before absorption. For many applications, it is certain that the ab-sorption occurs but this may take a very long time. Then, the quasi-stationary distribution gives an excellent measure of the long-term behavior of the system. Due to the non-linear structure of the quasi-stationary equations, it is usually impossible to obtain explicit expressions for the quasi-stationary distribution. However, there exists a number of helpful results including recursive methods, approximations and asymptotic analysis. The RE-distribution gives an alternative to measure the system dynamics before absorption, despite of how long the absorption time is. Since the RE-distribution is governed by linear equations, it can typically be evaluated more simply. The main problem of the RE-distribution concerns the practical difficulties for managing information about the initial distribution.

Several scenarios (classical SIS model, a linear growth model, classical SIR model, a finite SIR model with demography) are considered along the paper to investigate a two-fold objective: i) to consider the quasi-stationarity and the ratio of expectations as two conceptually different approaches for measuring the behavior of a biological system before reaching the absorbing states, and ii) to evaluate the possibility of using the RE-distribution as an approximation to the quasi-stationary distribution provided that the quasi-stationary regime has already been reached.

The study can be continued in several directions. For example, it would be interesting to get asymptotic expansions for the expectations E i [T j ] . If we allow absorption in one transition we arrive to the consideration of biological models with killing and catastrophes [START_REF] Coolen-Schrijner | Quasi-stationary distributions for a class of discrete-time Markov chains[END_REF][START_REF] Artalejo | Evaluating growth measures in an immigration process subject to binomial and geometric catastrophes[END_REF]. In the present paper, the applications of the RE-distribution were oriented to epidemic models but a forthcoming study might explore other important stochastic biological models including competition and predation [START_REF] Allen | An Introduction to Stochastic Processes with Applications to Biology[END_REF].

We also notice that the approach can be easily extended to the discretetime case. Similar arguments to those given in Section 2 lead to a system of equations for the number of visits to state j before absorption. In fact, the discrete counterpart follows by replacing the infinitesimal rates q ij by the onestep transition probabilities. Then, the analogue to equation (2.8) takes the form (I -P S T ) m j = e j , where P S T denotes the one-step transition probability sub-matrix corresponding to the transient states.

Fig. 3 .

 3 Fig. 3. States and transitions of the SIR model with demography

Table 1 .

 1 Comparing the characteristics of the SIS model with N = 100

		9284	5.0459	7.0314	19.7010	30.8972	48.93045166
		1.3152	3.9497	5.2233	8.9718	8.7257	7.23418980
		0.5141	0.1732	0.1091	0.0092	4.2168 × 10 -4	7.8067 × 10 -9
		1	1	1	20	32	49
	E[p (1) ] σ(p (1) )	1.9444 1.3312	5.5674 4.2725	8.1900 5.7378	23.2472 8.5418	33.3362 8.1589	50.00000001 7.07106775
	p (1) 1 p (1) m	0.5097 1	0.1501 1	0.0819 1	0.0022 24	4.4058 × 10 -5 34	3.2606 × 10 -10 50
	E[P 1 ]	1.4251	2.8837	4.0986	17.9807	30.7926	48.93044999
	σ(P 1 )	0.8602	2.7819	4.0950	9.7311	8.8501	7.23419422
	P 1 (1)	0.7264	0.4328	0.3356	0.0398	0.0013	1.5954 × 10 -8
	P m 1	1	1	1	20	32	49
	E[P N ]	17.7858	16.2558	16.0297	20.8561	30.9337	48.93045195
	σ(P N ) P N (1)	22.1172 0.1292	18.3569 0.0703	16.7672 0.0538	10.4969 0.0083	8.7764 4.1967 × 10 -4	7.23419022 7.8065 × 10 -9
	P m N	1	1	1	21	32	49
	E[ P]	1.8270	4.3226	6.4074	19.7563	30.8989	48.93045167
	σ( P)	1.0522	3.2331	4.6079	8.8935	8.7229	7.23418979
	P 1	0.4678	0.1736	0.1038	0.0087	4.2016 × 10 -4	7.8065 × 10 -9
	P m	1	1	1	20	32	49
	E[ P]	1.7397	4.2893	6.2809	19.6495	30.8970	48.93045166
	σ( P)	1.1181	3.5146	4.9186	8.9897	8.7258	7.23418980
	P 1	0.5672	0.2232	0.1383	0.0096	4.2217 × 10 -4	7.8065 × 10 -9
	P m	1	1	1	20	32	49
	E 1 [T ]	1.3765	2.3101	2.9789	25.1182	737.0896	62678430.2
	E N [T ]	7.7396	14.2113	18.5840	120.4356	2382.7997	128096977.8

Table 3 .

 3 Characteristics of the quasi-stationary distribution

		1.782729	3.358695	4.310522	5.197538	6.548307
	1	89.880561 1.157914	75.847244 2.531369	68.015593 3.200962	60.676247 3.745675	48.776272 4.447316
		(1, 92)	(1, 82)	(1, 76)	(1, 69)	(2, 55)
		1.868709	4.580918	7.014187	9.664507	13.639154
	λ = 0.5	95.195758 1.247307	85.992502 3.416681	78.088749 4.718968	69.072884 5.681166	53.392711 6.534392
		(1, 97)	(1, 94)	(1, 92)	(7, 74)	(14, 50)
		1.893701	5.279383	9.141594	13.925720	21.204501
	λ = 1.0	96.527713 1.274471	89.140401 3.913796	80.996214 5.704721	70.534363 6.792541	52.768775 7.085032
		(1, 98)	(1, 97)	(1, 97)	(13, 72)	(22, 51)
		1.920373	6.495080	14.149327	25.043827	39.970584
	λ = 5.0	97.753647 1.305275	92.160278 4.813970	82.903763 7.704323	69.670016 8.319430	51.334353 7.237475
		(1, 99)	(1, 99)	(1, 99)	(25, 70)	(40, 52)
		1.924342	6.744306	15.386401	27.732683	44.159848
	λ = 10.0	97.912692 1.310106	92.561982 5.003532	83.026944 8.137786	69.391751 8.526345	51.195365 7.238408
		(1, 99)	(1, 99)	(1, 99)	(27, 71)	(45, 51)

Table 4 .

 4 Characteristics of the P distribution

		1.761919	3.069476	4.085833	5.080829	6.600214
		90.365396	77.185894	69.161475	61.606620	48.957257
		0.972266	2.149048	2.854446	3.446865	4.186577
		(2, 90)	(3, 76)	(4, 68)	(5.61)	(7, 49)
		1.802639	4.563704	7.013047	9.768168	13.674764
	λ = 0.5	95.409072 1.021293	86.440183 3.140335	78.341271 4.523125	68.935081 5.568556	53.309221 6.477212
		(2, 95)	(5, 86)	(7, 78)	(10, 69)	(14, 53)
		1.815385	5.001877	9.140858	13.982033	21.206802
	λ = 1.0	96.899195 1.036569	89.731851 3.514131	81.066132 5.529311	70.447300 6.734425	52.764619 7.080823
		(2, 97)	(1, 97)	(9, 81)	(14, 71)	(22, 51)
		1.824432	6.076221	14.174662	25.049965	39.970584
	λ = 5.0	97.895901 1.048449	92.652233 4.294910	82.874674 7.556473	69.662724 8.310091	51.334353 7.237474
		(1, 99)	(1, 99)	(1, 99)	(25, 70)	(40, 52)
		1.825670	6.547522	15.398180	27.736629	44.159848
	λ = 10.0	98.026322 1.050229	92.785416 4.542308	83.013491 7.996709	69.387390 8.520852	51.195365 7.238408
		(1, 99)	(1, 99)	(13, 86)	(27, 71)	(45, 51)

Acknowledgements

This research was supported by the Spanish Government (Department of Education and Innovation) and the European Commission through project MTM 2008-01121. The authors would like to thank the referees for their helpful com-24 ments and suggestions which greatly improved the paper.