Marc-Michel Rohé

Nicolas Duchateau

Maxime Sermesant

Xavier Pennec

Combination of Polyaffine Transformations and Supervised Learning for the Automatic Diagnosis of LV Infarct

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Myocardial infarction occurs when blood flow to the heart muscle is lowered and the myocardial cells in the territory start dying. The local contractility is reduced and can lead, if prolonged, to severe remodelling of the heart to maintain physiological constraints [START_REF] Konstam | Left ventricular remodeling in heart failure. current concepts in clinical significance and assessment[END_REF]. The function of the heart is then impaired [START_REF] Bijnens | Investigating cardiac function using motion and deformation analysis in the setting of coronary artery disease[END_REF], and is no longer able to pump as efficiently as it used to, which might cause complications. Acute complications may include heart failure if the damaged heart is no longer able to pump blood adequately around the body. Therefore, a quantitative understanding of this pathology and how the heart function changes with an infarct is highly desired. Several methods for computer-aided diagnosis of infarct have already been developed using echocardiographic images of the heart coupled with pattern recognition algorithms [START_REF] Sudarshan | Automated Identification of Infarcted Myocardium Tissue Characterisation Using Ultrasound Images: A Review[END_REF] although none of the features used are explicitly related to physiological characteristics of cardiac function.

In this article, our goal is to classify between control subjects and patients with infarct in an automatic way, based on the STACOM statistical shape modeling labeled data [START_REF] Fonseca | The cardiac atlas project-an imaging database for computational modeling and statistical atlases of the heart[END_REF] which consist of a segmentation of the myocardium (both epi and endo) wall at end-diastole and at end-systole. These two categories of subjects may differ both in the shape of the heart and in the deformation along the cycle. Indeed, after an infarct the damaged region will tend to shrink and the deformation along the cycle will be lower. Similar studies have already been done with the same dataset as [START_REF] Zhang | Orthogonal shape modes describing clinical indices of remodeling[END_REF], which focuses on the shape differences between both population whereas we use both shape and motion features. Due to the complexity and high-dimensionality of these data, we try to quantify both shape and motion using a limited number of parameters, which we combine and use to compare patients and learn the main modes characterizing both populations.

The features of interest characterizing the shape of the patients consist of the regional thickness at both end-diastole and end-systole. We also use features representing the deformation along the cycle. Our approach relies on statistics on the motion of the heart between end-diastole and end-systole. We project the motion on the subspace of polyaffine transformations [START_REF] Arsigny | A fast and log-Euclidean polyaffine framework for locally linear registration[END_REF]. With these transformations, we can express a deformation with a limited number of parameters [START_REF] Mcleod | Spatio-temporal tensor decomposition of a polyaffine motion model for a better analysis of pathological left ventricular dynamics[END_REF]. We develop further the methodology by reducing the transformations to keep only the most relevant parameters.

Then, we test classical machine learning algorithms on our set of combined shape/motion parameters and compare the performance of each algorithm using cross-validation techniques. Validating the method with 10-fold cross-validation, we get results of 95% correct labeling on yet-unseen cases. In addition, our method notably highlights the relative importance of the different features for the classification of this population.

Extraction of features of interest through shape and motion dimensionality reduction

In this section, we introduce the first dimensionality reduction that is applied to the studied data (made of one segmentation at end-diastole and one at endsystole, each comprised of 1089 points both for the endocardium and the epicardium). It consists in a non-learning approach to project the data of these segmentations to a limited number of regional parameters representing motion and shape.

Polyaffine projection

Due to point to point correspondence of the meshes and prior registration, we already have an estimate of the displacement field φ mapping each point at end-diastole to the corresponding point at end-systole. Instead of looking at displacements fields, we choose to represent the cardiac motion by the stationary velocity fields (SVF) v such that v = log φ. Working with SVF allows to perform vectorial statistics on diffeomorphisms, while preserving the invertibility constraint, contrary to the Euclidian statistics on displacement fields.

In [START_REF] Arsigny | A fast and log-Euclidean polyaffine framework for locally linear registration[END_REF], the authors introduce the space of Log-Euclidean Polyaffine Transformations (LEPT). By defining K regions and smooth weights ω k (x), these transformations have the properties to describe locally affine deformations using few parameters while still being invertible. The polyaffine transformation is the weighted sum of these locally-affine transformations M k :

v poly (x) = K k=1 ω k (x)M k x.
In the case of cardiac motion, we have a standardized regional decomposition into the standard American Heart Association (AHA) 17 regions for the left ventricle. We define the weights ω k as normalized Gaussian functions around the barycenter xk of each region such that:

ωk (x) = exp κ 2 (x -xk ) T φ -1 k (x -xk ) , ω k (x) = ωk (x) N j=1 ωj (x)
.

If we gather the parameters of the polyaffine transformation into a large vector m such that m = vect(M 1 , ..., M K ). The parameters of the optimal projection of a Stationary Velocity Fields v onto the space of polyaffine transformations has an analytical solution [START_REF] Mcleod | Spatio-temporal tensor decomposition of a polyaffine motion model for a better analysis of pathological left ventricular dynamics[END_REF] m = m = Σ -1 b, which minimizes in the least-squares sense: In order to get interpretable parameters for each region, we choose to express them in a local coordinate system adapted to the geometry of the heart. If we call R = (O, e 1 , e 2 , e 3 ) the original Cartesian coordinate system, we define the local coordinate of the region k as

C(M 1 , ..., M K ) = Ω v poly (x) -v(x) 2 dx 1 2 (m -m) T Σ(m -m) -1 2 mΣ m.
R i = (O k , e k 1 , e k 2 , e k 3 )
where O k is the barycenter of the region (the red point in the enclosed figure), e 1 the radial vector (green vector), e 2 the longitudinal vector (purple vector) and e 3 the circumferential vector (blue vector). We can express the polyaffine parameters M = (R, T), where R is the 3 × 3 matrix of the rotational parameters and T is the translation, in this new frame through the equations:

R k = P -1 k R k P k T k = P -1 k (R k O k + T k ),
where P k is the transfer matrix from the base (e 1 , e 2 , e 3 ) to the base (e k 1 , e k 2 , e k 3 ). Then, the new expression of the parameters in this local coordinates system:

M =    s r a 1,2 a 1,3 t r a 2,1 s l a 2,3 t l a 3,1 a 3,2 s c t c    ,
can be related to physiological deformation. The 3 translation parameters correspond to the motion along the 3 local axes (radial, longitudinal, and circumferential) whereas the diagonal coefficients correspond to the strain along these directions.

We propose a method to further reduce the model by keeping only the 3 parameters of the motion and the 3 parameters of the strain. This defines a polyaffine projection that, when expressed in the local basis previously defined, has only these parameters not equal to zero. We first introduce the projection matrix Q which is a 12K × 6K matrix giving the relation between the 6K translation and diagonal parameters expressed in the local coordinates m L and the 12K parameters expressed in the original coordinates m, such that Qm L = m. When expressing m this way, we constrain it to be within the subspace spanned by Q. This subspace corresponds exactly to the polyaffine transformation whose non-diagonal and non-translation parameters are equal to zero in the local coordinates. The least-square minimization can now be rewritten as:

C(m) 1 2 (Qm L -m) T Σ(Qm L -m) - 1 2 mΣ m ∂C ∂m L = Q T Σ(Qm L -m) = 0 =⇒ m = Qm L = Q(Q T ΣQ) -1 Q T Σ m
For each of the 200 training data we compute the LEPT projection of the deformation field. We are able to parametrize the 3D displacement fields (made of 6534 parameters: 3 parameters for each of the 2178 points of the mesh) by only 6K = 102 polyaffine parameters. Despite this large reduction of dimensionality, these parameters explain on average more than 70% of the original displacement. Box-plots of each of the 6 parameters are shown in Figure 1, where the most discriminant parameters (p value < 0.001) are highlighted in bold. The radial displacement as well as the strain are significantly lower (in absolute value) for the infarcted subjects, which is consistent with what would be clinically expected. Similar differences can be seen for the longitudinal parameters. On the other side, the circumferential motion is less significant, mostly due to the fact that it is very hard to track it accurately with clinical images and therefore not reflected in the provided meshes.

Thickness parameters

On top of the polyaffine parameters that characterize the deformation of the heart during a cardiac cycle, we also introduce parameters representing the overall shape of the heart. We choose to study the thickness of the wall within each of the AHA zones at ED and ES. These parameters correspond to the initial and final stages of the transformation from ED to ES, and therefore complement the above-described parameters. We define the thickness as the local distance between endocardial points and their corresponding epicardial locations. These values are also averaged per AHA zone, and summarized in Fig. [START_REF] Bijnens | Investigating cardiac function using motion and deformation analysis in the setting of coronary artery disease[END_REF]. Significant differences are observed in the thickness of the myocardium wall at end-systole in most of the regions, especially near the apex, for the diseased patients with respect to the control group. On the other side, thickness ED diastole is less discriminant between both groups. Other parameters related to shape were considered (such as the height of the heart at ED/ES and the diameter at the base at ED/ES) but no significant differences between both populations were seen and therefore we do not use them for classification.

Dimensionality reduction of the parameters and classification

In this section, we use both polyaffine and thickness parameters previously introduced in order to classify between healthy and infarcted subjects. We use the machine-learning toolbox Scikit-Learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] to test a collection of standard stateof-the-art algorithms on our dataset and compare their performance in predicting yet unseen data. The features that serve to feed the tested learned algorithms were considered in four different ways: either polyaffine or thickness parameters separately (sets # 1 and # 2), or concatenated without normalization (# 3) or with normalization so that they have a mean of 0 and a variance of 1.

Learnt Dimensionality Reduction

Complementary to the a-priori reduction of dimensionality imposed by the polyaffine model and the use of 17 AHA regions, we also evaluated the influence of a second dimensionality reduction of the data both with a Principal Component Analysis (PCA) and a Principal Least Square (PLS) decomposition [START_REF] Rosipal | Overview and recent advances in partial least squares[END_REF] prior to the tested algorithm. PCA is designed to spread the data according to the main modes of variability and is known to be a useful dimension reduction pre-processing to prevent over-fitting and improve the performance of some machine-learning algorithm. PLS looks at modes of the input variables that correlate the most with an output variable (in our case the pathology label 0 or 1). Therefore, in contrast with PCA, the modes also correlate with our classification. In particular, Fig 3 summarizes the loadings of the first mode of the PLS with respect to each parameter. Notably, this can be used to assess which of the parameters is the most important for the classification. The radial parameters are the most prevalent, whereas both the circumferential parameters and the thickness at ED provide very little contribution to the first mode and therefore the classification.

Classification

All algorithms were tested with 10-fold cross validation on the dataset made of 200 patients. Fig. 4 summarizes the results of the different algorithms. Combining both sets of parameters improves the performance of most of the algorithms showing that these sets give different kind of information about the data. We also see that PLS regression, by preprocessing the data and orienting the modes of the input variables upto the best correlation with the pathology labels, improves the performance of all machine learning algorithms especially for Decision Tree and Nearest Neighbors. With more than 95% of correct labeling, SVM-SVC algorithm used on the PLS reduction with 5 modes is the method that performs the best. It is interesting to see that increasing the number of PLS modes further does not improve the classification. Our interpretation is that the subsequent modes of the PLS are not correlated to the classification and can therefore induce over-fitting of the data. We also tested the method with different cross-validation such as leave-one-out, 2-fold or 5-fold in order to see the robustness of the method with respect to the size of the training set and got similar performance.

Conclusion

In this paper, we evaluated the contribution of prior reduction of dimensionality to the classification of high-dimensional motion data. One of the assets of our work is an innovative methodology to project a motion on a reduced number of polyaffine parameters. We apply the methodology to classify a population and detect an infarct based on the segmentations at end-systole and end-diastole. Following the first dimensionality reduction given by the polyaffine parameters, we use traditional statistical reductions on our sets of parameters with PCA and PLS. Using 10-fold cross validation, we show that the resulting parameters have good predictive power with more than 95 % correct classification on 200 infarcted/control cases. We are also able to quantify the importance of each of the parameters in the classification. Notably, this provides insights into what is the main impact of an infarct both in terms of motion and shape.

4 )Fig. 1 .

 41 Fig. 1. Parameters of the polyaffine projection both for infarcted patients (red) and control subjects (blue). (Top row): radial parameters for both the diagonal parameters -representing strain -and the translation parameters -representing motion. (Middle row): longitudinal parameters. (Bottom row): circumferential parameters. In bold the most significant parameters (p-value < 0.001).

Fig. 2 .

 2 Fig. 2. Box-plots of the thickness of the myocardium wall per AHA region. In blue, the control set and in red the patients with infarcts. (Left): end-diastole. (Right): endsystole. In bold the most significant parameters (p-value < 0.001).

  Loadings of the first PLS mode showing the contribution of each of the parameters and each of the AHA zone. In green the most important parameters and in red the less important.

					STRAIN		TRANSLATION	THICKNESS
										END	END	SUM OVER ALL
				RADIAL LONG CIRC RADIAL LONG CIRC	DIASTOLE	SYSTOLE	PARAMETERS
			1	0.9%	0.3%	0.0%	2.1%	0.9%	0.5%	0.0%	0.2%	5%
			2	0.4%	0.9%	0.0%	2.3%	0.9%	0.1%	0.2%	0.2%	5%
		BASE	3 4	0.0% 0.0%	0.6% 0.0%	0.0% 0.2%	2.0% 1.5%	1.2% 1.3%	0.0% 0.0%	0.3% 0.0%	0.2% 0.0%	4% 3%
			5	0.1%	0.1%	0.3%	1.2%	1.2%	0.2%	0.0%	0.2%	3%
			6	0.5%	0.0%	0.0%	1.5%	1.0%	0.5%	0.0%	0.2%	4%
	AHA REGIONS	MIDDLE	7 8 9 10 11 12	1.9% 1.3% 0.9% 0.8% 0.9% 1.5%	1.4% 1.6% 1.6% 1.3% 0.5% 0.6%	0.2% 0.2% 0.1% 0.0% 0.0% 0.0%	2.3% 2.6% 2.4% 1.9% 1.5% 1.7%	0.8% 0.9% 1.1% 1.2% 1.2% 1.0%	0.5% 0.0% 0.1% 0.0% 0.0% 0.4%	0.0% 0.0% 0.0% 0.0% 0.0% 0.0%	0.8% 0.8% 0.6% 0.7% 0.7% 0.6%	8% 7% 7% 6% 5% 6%
		APICAL	13 14 15 16	1.9% 1.6% 1.4% 1.7%	1.3% 1.2% 0.7% 0.7%	0.4% 0.2% 0.0% 0.1%	2.5% 2.2% 1.6% 1.8%	0.8% 1.0% 1.1% 0.9%	0.1% 0.1% 0.0% 0.4%	0.3% 0.2% 0.1% 0.3%	1.6% 1.2% 1.1% 1.1%	9% 8% 6% 7%
		APEX	17	1.2%	1.1%	0.4%	1.5%	1.0%	0.0%	0.6%	1.1%	7%
	SUM OVER ALL REGIONS	17.0% 13.8% 2.2% 32.7% 17.8% 3.0%	2.1%	11.4%
	Fig. 3.

  Cross-validation results (10-fold) of the classification with respect to different state-of-the-art machine learning algorithms and different sets of input data. Combination of algorithms and parameters that have the best performance are shown in green whereas the worst are shown in red.

			Decision Tree	Random Forest	Logistic Regression	Nearest Neighbors	SVM Linear	SVM svc
		PolyAffine (PA)	79%	88%	85%	81%	89%	86%
		Thickness (TH)	76%	84%	86%	89%	85%	86%
		PA + TH	79%	91%	89%	90%	90%	92%
		Norm. PA + TH	77%	88%	92%	85%	90%	94%
		2 Modes	84%	87%	90%	89%	87%	89%
	PCA	5 Modes 10 Modes	84% 84%	87% 89%	90% 90%	87% 87%	88% 88%	93% 93%
		All modes	84%	80%	93%	87%	92%	96%
		2 Modes	91%	93%	94%	94%	93%	94%
	PLS	5 Modes 10 Modes	91% 91%	96% 97%	96% 96%	97% 94%	95% 95%	97% 96%
		All modes	91%	94%	94%	94%	93%	94%
	Fig. 4.

Ackowledgements The authors acknowledge the partial funding by the EU FP7funded project MD-Paedigree (Grant Agreement 600932)


Samuel Soubeyrand
email: samuel.soubeyrand@avignon.inra.fr

Construction of semi-Markov genetic-space-time SEIR models and inference Titre: Construction de modèles stochastiques génético-spatio-temporels et inférence

Keywords: Bayesian estimation, Genomic data, Spatiotemporal data, State-space model, Susceptible-Exposed-Infectious-Removed model Mots-clés : Estimation bayésienne, Données génomiques, Données spatiotemporelles, Modèle à espace d'état, Modèle Susceptible-Exposé-Infectieux-Retiré AMS 2000 subject classifications: 62H11, 62P10, 97K80

published or not. The documents may come

Construction of semi-Markov genetic-space-time SEIR models and inference

Introduction

Fast-evolving pathogens such as RNA viruses can cause human, animal and plant epidemics of high impact in developing and developed countries alike. For instance, hepatitis E causes 3 million acute cases each year worldwide (>60% in South and East Asia), of which 57,000 are fatal (World Health Organization, 2013). The 2014 Ebola outbreak, for which 4,507 cases were reported between 2013-12-30 and 2014-09-14, presented a real risk of continued expansion in fall 2014 (World Health Organization, 2014). During the foot-and-mouth outbreak in Great Britain in 2001, some 6 million animals were culled and the total cost was estimated to be in excess of £8 billion [START_REF] Anderson | Transmission dynamics and epidemiology of BSE in British cattle[END_REF][START_REF] Haydon | The UK foot-and-mouth disease outbreak -the aftermath[END_REF]. At the global scale and over three decades, the overall cost of sharka (infecting trees of the genus Prunus) was estimated to exceed C10 billion [START_REF] Cambra | Plum pox virus and the estimated costs associated with sharka disease[END_REF].

In order to minimize these social, environmental and economic costs, we need to most effectively control infectious diseases and thus to better understand how pathogens spread within host populations, yet this is something we know remarkably little about. Cases close together in their locations and timings are often thought to be linked, but timings and locations alone are usually consistent with many different scenarios of who infected whom. For fast-evolving pathogens such as the RNA viruses given as examples in the previous paragraph, pathogen genomic data can advantageously complete spatial and temporal data. The genome of such pathogens evolves so quickly relative to the rate that they are transmitted, that even over single short epidemics we can identify which hosts contain pathogens that are most closely related to each other. This information is valuable because when combined with the spatial and timing data it should help us infer more reliably who-transmitted-to-who over the course of a disease outbreak. However, doing this so that spatial, temporal and genetic data are appropriately combined remains a major statistical challenge. Furthermore, because sequencing genetic material has become so affordable, new statistical methods combining spatial, temporal and genetic data and estimating transmission trees will become very important for future epidemiology.

In recent years, the methodological challenge of reconstructing the dynamics of epidemics using (spatio-)temporal information and pathogen genetic information has been mainly addressed with two different but complementary approaches. The first approach is based on phylogeny / phylogeography, birth-death processes (often approximated by coalescent models) and BEAST [START_REF] Hall | Epidemic reconstruction in a phylogenetics framework: transmission trees as partitions[END_REF][START_REF] Lemey | Bayesian phylogeography finds its roots[END_REF][START_REF] Lemey | Phylogeography takes a relaxed random walk in continuous space and time[END_REF][START_REF] Pybus | Unifying the spatial epidemiology and molecular evolution of emerging epidemics[END_REF][START_REF] Rambaut | The genomic and epidemiological dynamics of human influenza A virus[END_REF][START_REF] Rasmussen | Inference for nonlinear epidemiological models using genealogies and time series[END_REF][START_REF] Stadler | Uncovering epidemiological dynamics in heterogeneous host populations using phylogenetic methods[END_REF]. This approach is used to estimate parameters related to the pathogen demographic process: the rate of spatial spread of the pathogen, the rate of switching between various host types, and the rate of evolution over time. The second approach is based on stochastic, spatiotemporal and evolutionary SEIR (Susceptible, Exposed, Infectious, Removed) models [START_REF] Jombart | Bayesian reconstruction of disease outbreaks by combining epidemiologic and genomic data[END_REF][START_REF] Morelli | A Bayesian inference framework to reconstruct transmission trees using epidemiological and genetic data[END_REF][START_REF] Mollentze | A Bayesian approach for inferring the dynamics of partially observed endemic infectious diseases from space-timegenetic data[END_REF][START_REF] Ypma | Unravelling transmission trees of infectious diseases by combining genetic and epidemiological data[END_REF][START_REF] Ypma | Genetic data provide evidence for wind-mediated transmission of highly pathogenic avian influenza[END_REF]. This modeling and inference approach combines heterogeneous and multi-scale processes and data: it links the epidemiological scale -or host population(s) scale-and the micro-evolutionary scale -or pathogen genome scale. The space-time-genetic SEIR models explicitly recognize the host population structure, and the epidemiological processes that govern the interaction of host and pathogen. They enable inferences to be made about epidemiological processes, specifically the transmission tree reflecting who infected whom (at the host, premise or population resolution), and other parameters (e.g. incubation durations, heterogeneity among hosts in susceptibility and force of infection, rates of evolution and spatial spread per transmission event, and sizes of infected populations).

In this article, we are interested in the second approach and the associated estimation algorithms allowing inference about transmission dynamics based on spatial, temporal and genetic data.

So far, the space-time-genetic SEIR approach has been developed in only a few articles (in generalist journals and journals of epidemiology and computational biology) that demonstrated that the approach is promising for reconstructing transmission trees and for estimating epidemiological parameters. However, this approach can be improved in several directions. One of these directions, on which we focus on in this article, concerns the inference method. In this article, we present the space-time-genetic SEIR model introduced by [START_REF] Mollentze | A Bayesian approach for inferring the dynamics of partially observed endemic infectious diseases from space-timegenetic data[END_REF] as a semi-Markov model (Section 2), we propose a new MCMC algorithm to estimate the model parameters and the transmission tree (Section 3; the originality of this algorithm lies in the reconstruction of transmitted pathogen sequences), and we assess the performance of the new algorithm with a simulation study (Section 4). In the simulation study, we especially assess the performance of the algorithm when pathogen genomic data are available for all sampled cases or for only a fraction of these cases. This introductory article gives probabilists and statisticians an invitation for enriching the class of space-time-genetic SEIR models and developing efficient inference algorithms in order to carry out high-impact analyses of infectious diseases due to viruses.

Model

The genetic-space-time SEIR model, presented below in Subsection 2.6, is a combination of a semi-Markov epidemic model and a Markovian evolutionary model. The following subsections shows how these submodels and the resulting genetic-space-time SEIR model are constructed.

Discrete-state, continuous-time Markovian SEIR model

Here, we consider a classical SEIR model describing the temporal dynamics of numbers of susceptible, exposed, infectious and removed individuals in a population affected by a pathogen [START_REF] Britton | Introduction to statistical inference for infectious diseases[END_REF] propose in this special issue of the journal an overview of statistical inference applied to epidemic models of this kind). Time is viewed as a continuous variable. Let S(t), E(t), I(t) and R(t) in N respectively denote the numbers of susceptible, exposed, infectious and removed individuals at time t ≥ 0. The sum of these quantities is equal to the instantaneous total size of the population T(t) ∈ N, i.e. S(t) + E(t) + I(t) + R(t) = T(t) for any time t ≥ 0.

In general, many different events can cause a change in the population pattern (S(t), E(t), I(t), R(t)), for instance the birth and death of susceptibles, the infection of susceptibles, the death of exposed individuals, the end of exposed stage (coinciding with the beginning of the infectious stage), the death in infectious individuals, and the end of infectious stage (coinciding with the beginning of removed stage).

For the sake of simplicity, in this article, we consider only three possible events, namely infection, end of exposed stage and end of infectious stage. Corresponding transition rates are provided in Table 1. In this model, the risk of infection is a combination of a basic risk whose rate is α 0 S(t), and an endogenous risk whose rate α 1 S(t)I(t)/T(t) is proportional to the number of infectious individuals in the population of interest. The basic risk may correspond to exogenous pathogen sources. For instance, in the case of zoonoses (i.e. diseases that can be transmitted from animals to humans) if the population of interest is the set of humans, the basic risk may correspond to animals infecting humans.

Description

Event Rate Infection S → S -1 & E → E + 1 α 0 S + α 1 SI/T End of exposed stage E → E -1 & I → I + 1 β E End of infectious stage I → I -1 & R → R + 1 δ I 2.2.
S k (t) ∑ j =k w k j I j (t)/T j (t).
In the latter rate, the weight w k j is a measure of the intensity of contacts between individuals of districts k and j. By convention, the intensity of contacts between individuals residing in the same district is equal to one. The spatial, discrete-state, continuous-time Markovian SEIR model considered in this section is defined by events and rates provided in Table 2.

Description Event Rate Infection S k → S k -1 & E k → E k + 1 α 0 S k + α 1 S k I k /T k + α 1 S k ∑ j =k w k j I j /T j End of exposed stage E k → E k -1 & I k → I k + 1 β E k End of infectious stage I k → I k -1 & R k → R k + 1 δ I k

Particular case: individual-based version of the model

Now, let us consider a particular case of the previous model: assume that for all k ∈ {1, . . . , n}, the size T(t) ≡ 1. Thus, districts are replaced by single individuals, and the model becomes an individual-based model where the dynamics of the epidemics is modeled at the individual resolution. In addition, values of S k (t), E k (t), I k (t) and R k (t) are in {0, 1} and their sum is equal to one whatever t. By assuming that each individual k ∈ {1, . . . , n} is located at x k in the planar space R 2 , events and rates shown in Table 2 can be re-written as in Table 3. The location x k can be viewed as the central or main location of k. We can see in Table 3 that the local endogenous risk disappeared from the expression of the rate of infection since a susceptible individual cannot infect himself (another justification is that S k (t)I k (t) = 0, ∀t ≥ 0). Moreover, the rate corresponding to the distant endogenous risk is now written α 1 ∑ j =k w(x jx k )I j where w is a kernel whose value depends on the relative locations of individuals k and j. In ecology and epidemiology, w is called a dispersal kernel, and w(x jx k ) measures the intensity of contact between individuals k and j. A classical and flexible class of dispersal kernels is formed by the power-exponential kernels parametrized by α 2 = (α 2,1 , α 2,2 ) and satisfying, for all x ∈ R 2 :

w(x) = α 2,2 2π(α 2,1 ) 2 Γ 2 α 2,2 exp - ||x|| α 2,1 α 2,2 , ( 1
)
where ||x|| is the Euclidean distance between the origine of the planar space and x. Thus, the measure of the intensity of contact between individuals k and j decreases with the distance separating the central locations of k and j. Note that in (1), the constant before the exponential ensures that the integral of w over R 2 is equal to one.

: 1 → 0 & E k : 0 → 1 α 0 + α 1 ∑ j =k w(x j -x k )I j End of exposed stage E k : 1 → 0 & I k : 0 → 1 β End of infectious stage I k : 1 → 0 & R k : 0 → 1 δ

Semi-Markov extension of the individual-based model

In the Markovian model presented in the previous subsection, the times spent by individuals in the exposed and the infectious states are exponentially distributed. Depending on the context, this may be viewed as an unrealistic assumption. For instance, the exposed duration, that corresponds to a latency or incubation duration, is usually not exponentially distributed but has a distribution with a mode away from zero (e.g. see [START_REF] Hampson | Transmission dynamics and prospects for the elimination of canine rabies[END_REF]. Semi-Markov models [START_REF] Barbu | Semi-Markov chains and hidden semi-Markov models toward applications -Their Use in Reliability and DNA Analysis[END_REF] offer a framework to handle non-exponential durations in some of the possible states. Thus, in this subsection, we introduce a semi-Markov model where durations in the exposed and infectious stages are independently drawn under gamma distributions (see Table 4). The draws are also independent from the duration in the susceptible stage.

Remark concerning the durations in the susceptible and removed states: the durations between successive changes of the susceptible state in the models of Subsection 2.1-2.4 are not exponentially distributed in general because the corresponding transition rates randomly vary in time (due to variations in the infectious states; see Appendix A). In the individual based models of Subsection 2.3-2.4, the removed states are absorbant and the corresponding durations are therefore infinite.

: 1 → 0 & E k : 0 → 1 α 0 + α 1 ∑ j =k w(x j -x k )I j End of exposed stage E k : 1 → 0 & I k : 0 → 1 Γ(β 1 , β 2 ) End of infectious stage I k : 1 → 0 & R k : 0 → 1 Γ(δ 1 ,

Markovian evolutionary model for a pathogen sequence

Now, suppose that the pathogen under consideration is an RNA virus that can evolves with time.

More specifically, we suppose that mutations can occur at s sites of the viral sequence between the four possible nucleobases that are adenine (A), cytosine (C), guanine (G) and uracil (U). We assume that mutations in different sites are independent but mutation rates vary as functions of the current nucleobases and the substituting nucleobases as in the 3-parameters Kimura substitution model [START_REF] Kimura | Estimation of evolutionary distances between homologous nucleotide sequences[END_REF]. Here, mutation and substitution are synonyms. Thus, at the s sites of the sequence under mutation, the mutation processes follow independent, discrete-state, continuoustime Markovian models that are defined by events and rates provided in Table 5.

Based on this model, [START_REF] Kimura | Estimation of evolutionary distances between homologous nucleotide sequences[END_REF] considered two sequences evolving from a hidden common ancestral sequence, they built a system of ordinary differential equations governing the probabilities of obtaining in the two sequences at any time all the possible pairs of nucleotides (AA, AC, AG, AU, CG, ...), and they solved the system analytically to obtain the expressions of the probabilities as a function of the time of evolution (further technical details are provided in [START_REF] Takahata | A model of evolutionary base substitutions and its application with special reference to rapid change of pseudogenes[END_REF].

In our framework, the times of evolution since a common ancestor of two observed sequences are generally different because the sequences are observed at varying times. However, because the mutation rates do not depend on the direction of the mutation (e.g. the same rate applies for A→G and G→A), the direction of the evolution does not matter and, formally, only the evolutionary time lag separating the two sequences matters. Thus, the expressions provided by [START_REF] Kimura | Estimation of evolutionary distances between homologous nucleotide sequences[END_REF] can be used in our framework to compute the expected proportions of the numbers of transitions, type-1 transversions, type-2 transversions and unchanged nucleobases over an evolutionary time lag ∆ separating two sequences. These expected proportions are:

ρ = (ρ 1 , ρ 2 , ρ 3 , ρ 4 ) = 1 4 1 -e 1 -e 2 + e 3 , 1 -e 1 + e 2 -e 3 , 1 + e 1 -e 2 -e 3 , 1 + e 1 + e 2 + e 3 , (2)
where e 1 = exp{-2(µ 1 + µ 2 )∆}, e 2 = exp{-2(µ 1 + µ 3 )∆}, e 3 = exp{-2(µ 2 + µ 3 )∆}, and µ 1 , µ 2 and µ 3 are the genetic substitution rates per nucleotide per day, for transitions, type-1 transversions and type-2 transversions, respectively.

In addition, we make the following distributional assumption: the numbers of observed transitions, type-1 transversions, type-2 transversions and unchanged nucleobases over an evolutionary time lag ∆ separating two sequences are distributed from a multinomial distribution, say P µ,s (• | ∆), with size equal to the length s of the observed sequence fragment and with the vector of probabilities ρ given by Equation (2). Thus, for any nonnegative integers m 1 , m 2 , m 3 , m 4 whose sum is s,

P µ,s {(m 1 , m 2 , m 3 , m 4 ) | ∆} = (s!) × ρ m 1 1 ρ m 2 2 ρ m 3 3 ρ m 4 4 (m 1 !) × (m 2 !) × (m 3 !) × (m 4 !) .

Genetic-space-time SEIR model

The genetic-space-time SEIR model, whose structure is illustrated by Figure 1, is obtained by combining the semi-Markov SEIR model of Subsection 2.4 and the Markovian evolutionary model of Subsection 2.5. The two models are combined under the following list of assumptions.

(t) = 1 or I k (t) = 1, is denoted by S k (t)
and is a vector of letters A, C, G and U. -When a susceptible individual k is infected by an infectious individual j at time t, the sequence S j (t) is transmitted to k, i.e. S k (t) = S j (t). -For any individual, mutations of the virus sequence during the exposed and infectious stages are assumed to be independent from the epidemiological dynamics. When an individual is removed, the sequence is fixed. This is illustrated in Figure 2 where the length of the sequence under mutation is s = 2. -Virus mutations in different infected individuals are assumed to be conditionally independent given the virus sequences at the infection times.

Disease reservoir ACCACGUC...

Exogenous infection

S j = 1 / / E j =1 ACCACGUC... / / I j = 1 ACGACGUC / / Distance-dependent transmission R j =1 ACGAUGUC... S k = 1 / / E k =1 ACGCGUC... / / I k = 1 / / R k =1 ACGACGCG... FIGURE 1.
Diagram illustrating the combination of the semi-Markov SEIR model and the Markovian evolutionary model. Individual j is infected by the disease reservoir with virus sequence "ACCACGUC...". Then, j becomes infectious and when j infects k, the sequence in j has evolved (C at the 3rd base mutated to G). Finally, the sequences in j and k independently evolve.

Soumis au Journal de la Société Française de Statistique File: seir.tex, compiled with jsfds, version : 2009/12/09 date: September 29, 2015 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q E=1 E=1 E=1 q q S=1 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q

E=1

I=1 I=1 I=1 I=1 I=1 I=1 I=1 I=1 I=1 I=1 I=1 I=1 I=1 I=1 I=1 I=1 AA AC AG AU CA CC CG CU GA GC GG GU UA UC UG UU q q q
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q

R=1 R=1 R=1 R=1 R=1 R=1 R=1 R=1 R=1 R=1 R=1 R=1 R=1 R=1 R=1 R=1 AA AC AG AU CA CC CG CU GA GC GG GU UA UC UG UU FIGURE 2.
Possible transitions for an individual in the genetic-space-time SEIR model, where the virus sequence is of length s = 2. Transitions from S = 1 (susceptible) to E = 1 (exposed), from E = 1 to I = 1 (infectious) and from I = 1 to R = 1 (removed) are irreversible, whereas transitions corresponding to substitutions in nucleobases are reversible.

Estimation methods

Data structure

We consider the same type of data than in [START_REF] Mollentze | A Bayesian approach for inferring the dynamics of partially observed endemic infectious diseases from space-timegenetic data[END_REF] who analyzed transmissions of rabies within a South African population of dogs, jackals and livestock. Thus, -Data are collected in a spatio-temporal observation window included in the whole spatial domain and in the whole time frame covered by the disease dynamics (which is either an epidemic or an endemic dynamics); -Among all the infected cases in the spatio-temporal observation window, only a subset is observed (thus, sources of infection can be unobserved, and the unobserved sources can be inside or outside the observation window); -The removed state corresponds to the death of the individual; -Sampled individuals are observed when they die, that is to say at the transition from I = 1 to R = 1; -Virus sequences that are available correspond to the states of the sequences at the death times (usually only a fragment of the sequence is available, the same fragment for all inividuals); -The central locations x 1 , . . . , x n of sampled individuals are assumed to be the locations at the death; -The sequence S exo is assumed to be known (in real cases, S exo can be reconstructed by a phylogenetic analysis of available genetic sequences: S exo is typically the reconstructed sequence of the most recent common ancestor; see [START_REF] Mollentze | A Bayesian approach for inferring the dynamics of partially observed endemic infectious diseases from space-timegenetic data[END_REF]. Compared to the amount of variables in the model, data are particularly sparse. Indeed, looking at Figure 1, data consist of the sequence of the disease reservoir and locations, times and sequences collected when observed individuals are in state R = 1. It has to be noted that, usually, only a 2, infected individuals and pathogen sequences are observed in one of the state R = 1 whereas they previously evolved in a high-dimension space of states (37 states in Figure 2; 1 + 2 × 4 s states in general, where s is the length of the observed sequence fragment).

Posterior distribution

Following [START_REF] Morelli | A Bayesian inference framework to reconstruct transmission trees using epidemiological and genetic data[END_REF] and [START_REF] Mollentze | A Bayesian approach for inferring the dynamics of partially observed endemic infectious diseases from space-timegenetic data[END_REF], we consider the joint posterior distribution p(J, T in f , L, D, θ | data) of the transmission tree J, infection times T in f = (T in f 1 , . . . , T in f n ), exposed (or latency) durations L = (L 1 , . . . , L n ), infectious durations D = (D 1 , . . . , D n ), and parameters θ that contains infection and dispersal parameters α = (α 0 , α 1 , α 2,1 , α 2,2 ), latency parameters β = (β 1 , β 2 ), infectiousness parameters δ = (δ 1 , δ 2 ), mutation parameters µ = (µ 1 , µ 2 , µ 3 ) and the date t exo of the exogenous sequence S exo .

The transmission tree J is a function from {1, . . . , n} to {0, 1, . . . , n} that states who infected whom: an observed individual i is infected by a pathogen source j = J(i) that is either another observed individual j ∈ {1, . . . , n}, j = i, or the disease reservoir (exogenous source) denoted by 0.

Data are removal times T end = (T end 1 , . . . , T end n ), central locations X = (x 1 , . . . , x n ) and observed sequences S end = {S 1 (T end 1 ), . . . , S n (T end n )}. The posterior distribution is:

p(J, T in f ,L, D, θ | data) = p(J, T in f , L, D, θ | S end , T end , X, S exo ) ∝ p(S end | J, T in f , L, D, θ , T end , X, S exo )p(J, T in f , L, D, θ | T end , X, S exo ) = p(S end | J, T in f , L, D, θ , T end , X, S exo )p(J, T in f | L, D, θ , T end , X, S exo ) × p(L, D | θ , T end , X, S exo )p(θ ), (3)
where ∝ means "proportional to" (the multiplicative constant does not depend on the unknowns (J, T in f , L, D, θ )), p(S end | J, T in f , L, D, θ , T end , X, S exo ) is called the genetic likelihood, p(J, T in f | L, D, θ , T end , X, S exo ) is called the transmission likelihood, p(L, D | θ , T end , X, S exo ) is the distribution of latency and infectious durations and p(θ ) is the prior distribution of parameters.

The prior distribution of parameters will be specified in Subsection 4.2. Based on Subsection 2.4, the distribution of latency and infectious durations is simply the product of gamma probabilities: In Subsection 3.3, we detail the transmission likelihood and show how the incompleteness of epidemiological data, i.e. the missing infecting hosts, is handled. In Subsections 3.4, 3.5 and 3.6, we detail the genetic likelihood and show how the incompleteness of genetic data, i.e. the missing pathogen sequences, is handled. The genetic likelihood can be formally written as a function of the unobserved virus sequences at the infection times S k (T in f k ), k = 1, . . . , n. To avoid 2014) replaced the genetic likelihood by a pseudo-likelihood in their MCMC algorithm; see Subsection 3.5. In Subsection 3.6, we propose a new approximation of the genetic likelihood whose performance will be assessed in Section 4.

p(L, D | θ , T end , X, S exo ) = p(L, D | θ ) = I ∏ i=1 γ(L i ; β 1 , β 2 )γ(D i ; δ 1 , δ 2 ),

Transmission likelihood

Following [START_REF] Mollentze | A Bayesian approach for inferring the dynamics of partially observed endemic infectious diseases from space-timegenetic data[END_REF], the transmission likelihood p(J, T in f | L, D, θ , T end , X, S exo ) can be written:

p(J, T in f |L, D, θ , T end , X, S exo ) = p J(1), T in f 1 | L, D, θ , T end , X × I ∏ i=2 p J(i), T in f i | J{1 : (i -1)}, T inf 1:(i-1) , L, D, θ , T end , X , (4)
where the index i is sorted with respect to increasing infection times T in f i , J{1 : (i -1)} = (J(1), . . . , J(i -1)) for i > 1, T inf 1:(i-1) = (T in f 1 , . . . , T inf i-1 ) for i > 1, and by assuming that the transmission dynamics does not depend on the exogenous sequence S exo .

Each host has the same chance (1/I) to be infected first (by an external source J(1) = 0), and its infection time is assumed to be less than or equal to the first observation time (min{T end }):

p J(1), T in f 1 | L, D, θ , T end , X = 1 I × 1(T in f 1 ≤ min{T end }),
where 1 is the indicator function. Subsequent infections (i.e. for i > 1) occur with the following probabilities derived from assumptions made in Subsections 2.4 and 2.6:

p J(i), T in f i | J{1 : (i -1)}, T inf 1:(i-1) , L, D, θ , T end , X = exp -α 0 (T in f i -T in f 1 ) - T in f i T in f 1 i-1 ∑ j=1 α 1 1(T in f j + L j ≤ t ≤ T end j )w(x j -x i )dt × α 0 1{J(i) = 0} + α 1 1(T inf J(i) + L J(i) ≤ T in f i ≤ T obs J(i) )w(x J(i) -x i )1{J(i) > 0}
where the exponential term is the probability that host i has not been infected between times T in f 1 and T in f i , and the second term is the probability density that host i has been infected by J(i) at time T in f i . Here, if J(i) > 0 the source is observed, while the source is external to the dataset (an introduction) if J(i) = 0. α 0 is the infection strength of the exogenous sources, assumed to be constant in time and space, α 1 is the infection strength of an observed source, and w is the parametric dispersal kernel specified in Equation (1).

Genetic likelihood

By assuming that the evolution of pathogen sequences does not depend on the transition from exposed to infectious states and does not depend on individual locations, the genetic likelihood p(S end | J, T in f , L, D, θ , T end , X, S exo ) reduces to p(S end | J, T in f , θ , T end , S exo ). This distribution was written by [START_REF] Mollentze | A Bayesian approach for inferring the dynamics of partially observed endemic infectious diseases from space-timegenetic data[END_REF] as a sum over all possible transmitted sequences S in f i = S i (T in f i ) at times T in f i for all i such that J(i) > 0:

p(S end | J, T in f , θ , T end , S exo ) = ∑ {S in f i ∈S: i=1,...,n and
J(i)>0} n ∏ i=1 P µ,s {M(S end i , S † prec(i,obs) ) | ∆ = T end i -T † prec(i,obs) } ×    n ∏ i=1 J(i)>0 P µ,s {M(S in f i , S * prec(i,in f ) ) | ∆ = T in f i -T * prec(i,in f ) }         . (5)
In Equation ( 5), the first series of factors accounts for the probabilities of the number of substitutions between an observed sequence and the immediately preceding unobserved, transmitted sequence or the sequence of the exogenous source (last double arrows on every horizontal lines of Figure 3A). The second series of factors accounts for the probabilities of the number of substitutions between each transmitted sequence and the transmitted sequence immediately preceding in time or the sequence of the exogenous source (every double arrows except the last one on every horizontal lines of Figure 3A).

Let us now detail terms in Equation ( 5): S is the set of all possible sequences (the size of S is 4 s , where s is the length of the sequence fragment); M(S , S) is the vector of the numbers of transitions, type-1 transversions, type-2 transversions and unchanged nucleobases between S and S ; P µ,s {M(S , S) | ∆ = T -T } is the probability given by the multinomial distribution in Subsection 2.5.

The subscript prec(i, obs) can take two types of values. First case: if J(i) = 0 and i did not infected any other observed host, then prec(i, obs) = exo, S † prec(i,obs) = S exo and T † prec(i,obs) = t exo . Second case: if J(i) > 0 or if i infected another observed host, then prec(i, obs) denotes the host whose node of infection belongs to the tree path from the root of the tree to the observation of i (at time T end i ) and whose infection is just preceding the observation of i (S † prec(i,obs) is the transmitted sequence S in f prec(i,obs) at the infection time T † prec(i,obs) = T in f prec(i,obs) ). The node of infection of a given host k is defined as the point on the tree at which "the branch leading to the observation of k" and "the branch leading to the observation of the infecting host J(k)" diverged. The tree path from one point of the tree to another is defined as the most direct broken line on the graph connecting the two points. In the second case, if J(i) > 0 and i did not infect any other host, then prec(i, obs) is i itself.

The subscript prec(i, in f ) can also take two types of values. First case: if J(J(i)) = 0 and J(i) did not infect any other observed host before the infection of i at T in f i , then prec(i, in f ) = exo, S * prec(i,in f ) = S exo and T * prec(i,in f ) = t exo . Second case: if J(J(i)) > 0 or if J(i) infected another observed host before the infection of i at T in f i , then prec(i, in f ) denotes the host whose node of infection belongs to the tree path from the root of the tree to the infection of i (at time T in f i ) and whose infection is just preceding the infection of i (S * prec(i,in f ) is the transmitted sequence S in f prec(i,in f ) at the infection time T * prec(i,in f ) = T in f prec(i,in f ) ).

Genetic pseudo-likelihood

To reduce the complexity of the inference algorithm, Mollentze et al. ( 2014) used a conditional pseudo-distribution of S end , noted p(S end | J, T in f , θ , T end , S exo ), instead of the exact conditional distribution given by Equation ( 5). The conditional pseudo-distribution does not depend on the extra latent vectors {S in f i : i = 1, . . . , n, J(i) > 0} appearing in Equation ( 5). With index i being sorted with respect to increasing infection times T in f i , the distribution p(S end | J, T in f , θ , T end , S exo ) can be written:

p(S end | J, T in f , θ , T end , S exo ) =p(S end 1 | J, T in f , θ , T end , S exo ) × n ∏ i=2 p(S end i | S end 1:(i-1) , J, T in f , θ , T end , S exo ), (6)
where S end 1:(i-1) is the set of observed sequences from hosts 1, . . . , i -1. For the first infected host,

p(S end 1 | J, T in f , θ , T end , S exo ) = P µ,s {M(S end 1 , S exo ) | ∆ = T end 1 -t exo }.
For the other hosts infected by the exogenous source (i.e. for i > 1 such that J(i) = 0),

p(S end i | S end 1:(i-1) , J, T in f , θ , T end , S exo ) = P µ,s {M(S end i , S exo ) | ∆ = T end i -t exo }.
For hosts infected by observed hosts (i.e. for i > 1 such that J(i) > 0), we replaced the conditional probability p(S end i | S end 1:(i-1) , J, T in f , θ , T end , S exo ) of S end i given sequences S end j ( j = 1, . . . , i -1) by the product, up to a power, of the conditional probabilities of S end i given each sequence S end j such that j ∈ 1, . . . , i -1 and j is in the transmission chain leading to i (the latter condition is mathematically written: ∃m ∈ N * , J m (i) = j, where J m consists of composing J with itself m times):

    i-1 ∏ j=1 ∃m∈N * ,J m (i)= j P m,s {M(S end i , S end j ) | ∆ =| T end i -T in f div(i, j) | + | T end j -T in f div(i, j) |}     1/η i , (7)
where η i is the number of terms in the product, T end div(i, j) denotes the infection time at which the chain of infection leading to i and the chain of infection leading to j diverged (T in f div(i, j) is one of the latent variables in T in f ) and

∆ =| T end i -T in f div(i, j) | + | T end j -T in f div(i, j)
| is the evolutionary duration separating the observation of S end i and S end j . The use of the power 1/η i is a way to get a quantity homogeneous to a single probability and not to a product of probabilities whatever the length of the transmission chain leading to i. Therefore, the hosts have similar weights in the pseudo-distribution given below. The computation of Equation ( 7) is illustrated by Figure 3B. Thus, the conditional pseudo-distribution of S end , or in other words the genetic pseudolikelihood, satisfies:

p(S end | J, T in f , θ , T end , S exo ) = n ∏ i=1 J(i)=0 P µ,s {M(S end i , S exo ) | ∆ = T end i -t exo } × n ∏ i=1 J(i)>0     i-1 ∏ j=1 ∃m∈N * ,J m (i)= j P m,s {M(S end i , S end j ) | ∆ =| T end i -T in f div(i, j) | + | T end j -T in f div(i, j) |}     1/η i . (8)

Genetic approximate likelihood

The genetic likelihood given by Equation ( 5) is a sum over all possible transmitted sequences S in f i for i such that J(i) > 0. In the genetic approximate likelihood, we use a fixed value Sin f i for S in f i that is deterministically reconstructed conditionally on J, S exo and S end . The reconstruction is based on the parsimony principle commonly used in phylogenetics: the most parsimonious reconstruction of {S in f i : i = 1, . . . , n, J(i) > 0} is the one that requires the fewest evolutionary changes (i.e. the fewest substitutions of nucleobases; see [START_REF] Tuffley | Links between maximum likelihood and maximum parsimony under a simple model of site substitution[END_REF], for a formal definition). The genetic approximate likelihood is given by: p(S end | J, T in f , θ , T end , S exo )

= n ∏ i=1 P µ,s {M(S end i , S † prec(i,obs) ) | ∆ = T end i -T † prec(i,obs) } ×    n ∏ i=1 J(i)>0 P µ,s {M( Sin f i , S * prec(i,in f ) ) | ∆ = T in f i -T * prec(i,in f ) }    , (9)
where S † prec(i,obs) and S * prec(i,in f ) are computed like S † prec(i,obs) and S * prec(i,in f ) in Subsection 3.4 except that S in f i is replaced by Sin f i for all i such that J(i) > 0. The set of sequences { Sin f i : i = 1, . . . , n, J(i) > 0} is computed by applying Algorithm 2 that calls Algorithm 1 several times.

Algorithm 1 computes one of the most parsimonious sequence for the node of bifurcation of an ancestral sequence G anc into two sequences G (1) and G (2) (there can be several most parsimonious sequences for the node; see Discussion section). In our setting (and in Algorithm 2), nodes of bifurcation coincide with infection events. The function that provides G node by carrying out calculations in Algorithm 1 is denoted by G :

(G anc , G (1) , G (2) ) → G(G anc , G (1) , G (2) ).
Algorithm 2 shows how the fixed value { Sin f i : i = 1, . . . , n, J(i) > 0} are obtained given J, S exo and S end . Note that the operation G anc ⇐ Sin f j that uses one of the output of the algorithm is feasible because index i is sorted with respect to increasing infection times T in f i (i.e. Sin f j has been already computed when the operation G anc ⇐ Sin f j is performed).

Remark: Algorithm 2 is not demonstrated to yield the most parsimonious reconstruction of the transmitted sequences. However, it is expected to yield a relatively parsimonious reconstruction because it is based on successive calls of Algorithm 1, which leads to the most parsimonious sequence for the node of bifurcation of an ancestral sequence G anc into two sequences G (1) and G (2) .

Algorithm 1 Most parsimonious reconstruction of the sequence G node = (g node 1 , . . . , g node s ) of the node of bifurcation from an ancestral sequence G anc = (g anc 1 , . . . , g anc s ) to two sequences

G (1) = (g (1) 1 , . . . , g (1) s ) and G (2) = (g (2) 1 , . . . , g (2)
s ). These sequences, that are shown on the diagram below on the right of the algorithm, are vectors of size s whose components take values in a finite countable set.

for b = 1, . . . , s do if g (1) b = g (2) b then g node b ⇐ g (1) b else g node b ⇐ g anc b end if end for G (1) G anc / / G node 4 4 i i i i i i i i i i i i i * * U U U U U U U U U U U U U G (2)
Algorithm 2 Reconstruction of sequences {S in f i : i = 1, . . . , n, J(i) > 0} conditionally on J, S exo and S end . Index i is sorted with respect to increasing infection times T in f i .

for j = 1, . . . , (n -1) do if J( j) = 0 then G anc ⇐ S exo else G anc ⇐ Sin f j end if for i = 1, . . . , (n -1) do if J(i) = j then G (1) ⇐ S end j G (2) ⇐ S end i Sin f i ⇐ G(G anc , G (1) , G (2) ) G anc ⇐ Sin f i end if

MCMC

To estimate parameters and latent variables, in particular the transmission tree J, [START_REF] Mollentze | A Bayesian approach for inferring the dynamics of partially observed endemic infectious diseases from space-timegenetic data[END_REF] built an MCMC algorithm that samples in:

p(J, T in f ,L, D, θ | data) ∝ p(S end | J, T in f , L, D, θ , T end , X, S exo )p(J, T in f | L, D, θ , T end , X, S exo ) × p(L, D | θ , T end , X, S exo )p(θ ), (10)
Equation ( 10) is analogue to Equation ( 3) that gives the posterior distribution of parameters and latent variables, except that the genetic likelihood (Eq. ( 5)) is replaced by the genetic pseudolikelihood (Eq. ( 8)).

Here, we propose to estimate parameters and latent variables with an MCMC algorithm that samples in:

p(J, T in f ,L, D, θ | data) ∝ p(S end | J, T in f , L, D, θ , T end , X, S exo )p(J, T in f | L, D, θ , T end , X, S exo ) × p(L, D | θ , T end , X, S exo )p(θ ), (11)
Equation ( 11) is analogue to Equation ( 3), except that the genetic likelihood (Eq. ( 5)) is replaced by the genetic approximate likelihood (Eq. ( 9)).

As explained above, sampling in p(J, T in f , L, D, θ | data) instead of sampling in the exact posterior avoids us to handle the latent transmitted viral sequences. In contrast, sampling in p(J, T in f , L, D, θ | data) requires us to handle these latent sequences but, using Algorithm 2, they are deterministically calculated given the values of J, S exo and S end . In both approaches, avoiding to stochastically update the transmitted viral sequences along the MCMC algorithm (with a Metropolis-Hastings sampler for example) significantly reduces the complexity of the algorithm. Note that computation costs are comparable when the pseudo-likelihood and the approximate likelihood are used: sampling in p(J, T in f , L, D, θ | data) only took us about 1.05 times the duration required for sampling in p(J, T in f , L, D, θ | data).

For the application presented in the next section, we used both MCMC algorithms (the one based the pseudo-likelihood and the one based on the approximate likelihood). Proposal distributions and initialization were approximately similar than those used by [START_REF] Mollentze | A Bayesian approach for inferring the dynamics of partially observed endemic infectious diseases from space-timegenetic data[END_REF]. It has to be noted that the transmission tree J was initialized such that J(i) = 0 for all i ∈ {1, . . . , n} (i.e. all individuals are initially infected by the exogenous source). Each MCMC algorithm was run for 100,000 iterations. Following a burn-in of 10,000 iterations, we sampled every 250th iterations.

Application

We carried out a simulation study to assess the efficiency of the two MCMC algorithms presented in Subsection 3.7. In this article, we assessed the algorithm efficiency by focusing on the estimation of the transmission tree J. We also evaluated how this efficiency is decreased when a fraction of genetic data are missing (i.e. when pathogen is not collected on some of the observed individuals and, therefore, not sequenced). Regarding the estimation of model parameters that is not studied here, [START_REF] Morelli | A Bayesian inference framework to reconstruct transmission trees using epidemiological and genetic data[END_REF] and [START_REF] Mollentze | A Bayesian approach for inferring the dynamics of partially observed endemic infectious diseases from space-timegenetic data[END_REF] provided extensive results on the coverage of parameters by posterior intervals provided by our approach.

Simulation model

For each simulation, the epidemic is initiated at time zero with one infected host localized at the origin (0, 0) and 119 susceptible hosts uniformly and randomly localized in the [0.0, 0.3] × [0.0, 0.1] rectangle. At time zero, the sequence of the virus in the infected host purely consists of A nucleobases. The incubation parameters are β 1 = 50 (mean) and β 2 = 5 (standard deviation). The infectious duration parameters are δ 1 = 10 (mean) and δ 2 = 1 (standard deviation). The unit of β 1 , β 2 , δ 1 and δ 2 is an arbitrary time unit that is typically the day. The infection strength of each infectious host is α 1 = 10 4 . The basic risk parameter is α 0 = 0. The dispersal parameters are α 2,1 = 0.002 (scale) and α 2,2 = 1 (shape). The substitution rates are µ 1 = µ 2 = µ = 2 × 10 -5 (in substitutions per nucleotide per time unit). This value is rather high but is included in current values observed for RNA viruses [START_REF] Hanada | A large variation in the rates of synonymous substitution for rna viruses and its relationship to a diversity of viral infection and transmission modes[END_REF][START_REF] Jenkins | Rates of molecular evolution in rna viruses: a quantitative phylogenetic analysis[END_REF]. Only a fraction of the n = 120 hosts is sampled: hosts with x-coordinate larger than 0.2 have a probability equal to 3/4 to be sampled, hosts with x-coordinate lower than or equal to 0.2 are not sampled (the sampling region is reduced to the square [0.2, 0.3] × [0.0, 0.1]). For genetic data, a sequence fragment of length s = 800 sites is observed. The values β 1 = 50, δ 1 = 10, µ 1 = µ 2 = µ 3 = 2 × 10 -5 and s = 800 lead, in expectation, to 0.96 substitutions per infected host over the observed sequence fragment. Figure 4 presents an example of simulated epidemic.

The simulation model is not exactly the same than the model used for the inference. Indeed, there is a difference concerning infection from external sources. In the simulation model, the external sources are handled in a realistic manner: the external sources are infectious hosts within or outside the sampling region generating infectious risks that are local in time and space. In the inference model, exogenous infections are due to a single source that has a constant infection strength (measured by α 0 ) in time and space and that is characterized by the sequence S exo at time t exo . This trick allows us to not have to explicitly deal with multiple unobserved sources (and therefore supplementary latent variables) in the estimation algorithm. It has to be noted that the parameter α 0 has not the same meaning in the simulation model and the inference model.

Prior distributions

Independent prior distributions were used for all parameters. Exponential priors with mean values 10 6 were chosen for α 0 and α 1 . An exponential prior with mean value 1 was chosen for α 2,1 , while an informative gamma prior distribution with mean 1 and standard deviation 1 was specified for α 2,2 . Doing so allows classical kernels to have a non-negligible weight, in particular the thin-tailed normal kernel (α 2,2 = 2), the exponential kernel (α 2,2 = 1) and the fat-tailed kernel corresponding to α 2,2 = 0.5.

For the parameters governing incubation and infectious periods, we used very narrow prior distributions centered around true parameter values. Thus, we consider a situation where accurate information about these quantities are available (e.g. see [START_REF] Hampson | Transmission dynamics and prospects for the elimination of canine rabies[END_REF], for such information in the case of canine rabies). Narrow prior distributions were used to prevent the inference algorithm from creating direct connections by extending the incubation and infectious periods when one or more intermediate cases have not been sampled (Note that in a situation where all infected cases are observed, the parameters governing incubation and infectious periods can be estimated with our approach [START_REF] Morelli | A Bayesian inference framework to reconstruct transmission trees using epidemiological and genetic data[END_REF]. A narrow gamma prior with mean 50 and standard deviation 0.01 was specified for β 1 . A narrow gamma prior with mean 5 days and standard deviation 0.01 was specified for β 2 . A narrow gamma prior with mean 10 and standard deviation 0.01 was specified for δ 1 . A narrow gamma prior with mean 1 and standard deviation 0.01 was specified for δ 2 .

Exponential prior distributions with mean parameter m = 2 × 10 -5 substitutions per nucleotide per day were used for the substitution rates µ 1 , µ 2 and µ 3 . A normal prior distribution with mean 0 and standard deviation 10 was specified for t exo .

Results

Fifty data sets were independently simulated using the procedure described above. For each data set, we applied four different estimation algorithms:

-in the first one, we used the genetic pseudo-likelihood and all genetic data; -in the second one, we used the genetic approximate likelihood and all genetic data; -in the third one, we used the genetic approximate likelihood and 50% of genetic data (50% of the available sequences were randomly selected and used for the inference, the other sequences were ignored); -in the fourth one, we used the genetic approximate likelihood and 25% of genetic data (25% of the available sequences were randomly selected and used for the inference, the other sequences were ignored). It has to be noted that the genetic sampling effort (which varies) is different from the epidemiological sampling effort (which is constant in expectation): the sampling scheme described in Subsection 4.1 leads to observe about 3/4 of the infected hosts located within the study region; then, the virus is sampled and sequenced for either 100%, 50% or 25% of the observed hosts.

Comparing the two first algorithms will allow us to assess which of the approaches based on the pseudo-likelihood and the approximate likelihood is more accurate for reconstructing the transmission tree in the simulation setting described in Subsection 4.1.

Comparing the three last algorithms will allow us to assess how the reconstruction of the transmission tree relies on the genetic data, by keeping constant the other components of the reconstruction (i.e. the epidemiological data and the inference algorithm).

The algorithms were compared by assessing their ability to correctly estimate endogenous and exogenous transmissions. Endogenous transmissions designate direct transmissions between observed cases whereas exogenous transmissions designate infections of observed cases by unobserved cases. Three criteria were used to compare the performance of the algorithms. They are defined as follows, for each simulation and each algorithm:

-Average posterior probability of true endogenous transmissions:

C 1 = 1 ∑ n i=1 1{J true (i) > 0} n ∑ i=1 1{J true (i) > 0} p{J(i) = J true (i) | data},
where J true is the true transmission tree that was simulated, p{J(i) = J true (i) | data} is the MCMC-based assessment of the posterior probability that the infecting source of i is the z) is the z-th state of J in the posterior sample obtained from the MCMC, and Z is the size of the posterior sample of (J, T in f , L, D, θ ); -Average posterior probability of true exogenous transmissions:

p{J(i) = J true (i) | data} = 1 Z Z ∑ z=1 1{J (z) (i) = J true (i)}}, J (
C 2 = 1 ∑ n i=1 1{J true (i) = 0} n ∑ i=1 1{J true (i) = 0} p{J(i) = J true (i) | data};
-Average posterior probability of false exogenous transmissions:

C 3 = 1 ∑ n i=1 1{J true (i) > 0} n ∑ i=1 1{J true (i) > 0} p{J(i) = 0 | data},
where p{J(i) = 0 | data} is the MCMC-based assessment of the posterior probability that the infecting source of i is an external source,

p{J(i) = 0 | data} = 1 Z Z ∑ z=1 1{J (z) (i) = 0}}.
Table 6 provides the means and standard deviations of C 1 , C 2 and C 3 computed from the 50 simulations and using the four estimation algorithms. Using the genetic approximate likelihood instead of the pseudo-likelihood allowed us to significantly improve the identification of true endogenous sources when all genetic data were used (the mean of C 1 was increased from 0.75 to 0.82; p-value of pairwise t-test: 0.018).

When the genetic sampling effort was decreased from 100% to 50% and 25%, the correct identification of both endogenous and exogenous sources was more uncertain but remained relatively high. Indeed, in average, there were 31 observed hosts in the simulated data sets and consequently, for each infected host there were 30 possible endogenous sources and 1 possible exogenous source. Therefore, obtaining mean values for C 1 and C 2 equal to 0.36 and 0.72 indicates that, in average, the true source is identified with a relatively high probability. We can notice that the decrease in the identification of true endogenous transmissions is mostly counterbalanced by an increase in the incorrect detection of exogenous transmissions. Thus, in general, when an endogenous source is not identified correctly, the algorithm tends to assign the infection to an external source but not to an incorrect endogenous source.

Discussion

In this article, we presented an extension of stochastic SEIR models that combines (i) an individualbased, spatial, semi-Markov SEIR model for the spatio-temporal dynamics of the pathogen, and (ii) a Markovian evolutionary model for the temporal evolution of genetic sequences of the pathogen. The resulting model is a state-space model including many latent variables. The high 6. Means and standard deviations of average posterior probabilities of true endogenous transmissions (C 1 ), true exogenous transmissions (C 2 ) and false exogenous transmissions (C 3 ). Means and standard deviations were computed from 50 simulations & estimations. Endogenous transmissions designate direct transmissions between observed cases whereas exogenous transmissions designate infections of observed cases by unobserved cases. Four estimation algorithms were applied, one based on the genetic pseudo-likelihood and three based on the genetic approximate likelihood with three different levels of genetic sampling effort. The genetic sampling effort is the percentage of observed cases for which the virus was sequenced and used in the inference.

Method

Pseudo We assessed this new algorithm with respect to its capacity to reconstruct transmission trees (i.e. who infected whom). This new algorithm outperforms the one proposed by [START_REF] Mollentze | A Bayesian approach for inferring the dynamics of partially observed endemic infectious diseases from space-timegenetic data[END_REF]. We also assessed the performance of the algorithm when pathogen sequences are observed for only a fraction of sampled hosts. The results that we obtained show that including more genetic data significantly improves the reconstruction of transmission trees. It has to be noted that similar studies could be carried out to investigate, for example, the influence of spatial data, contact tracing information, substitution rate heterogeneity and dissemination capacities in the reconstruction of transmission trees, but also in the estimation of key epidemic parameters and latent variables such as unobserved infection times. These influence analyses were out of the scope of the present introductory article but they are obviously of major interest. Several directions can be explored to improve the space-time-genetic SEIR approach. One of these directions concerns the computation / approximation of the genetic likelihood. We based our approximation on the reconstruction of transmitted genetic sequences based on the parsimony principle. While the parsimony principle has drawbacks in phylogenetic analyses [START_REF] Brocchieri | Phylogenetic inferences from molecular sequences: Review and critique[END_REF], it can be viewed as a nonparametric approach that is robust to model misspecifications [START_REF] Kolaczkowski | Performance of maximum parsimony and likelihood phylogenetics when evolution is heterogeneous[END_REF]. One of its main interests is its low computation cost with respect to the costs of probabilistic alternatives reviewed in [START_REF] Brocchieri | Phylogenetic inferences from molecular sequences: Review and critique[END_REF]. However, we can improve how the parsimony principle is applied in our approach. Let us give below an example of possible improvement (notations are those of Section 3.6). If g

b (resp. g

b ). Thus, the parsimony principle does not lead to a unique solution. The criterion using the differences in nucleobases could be augmented by a criterion on temporal lengths of branches joining G anc -G node , G node -G (1) , and G node -G (2) . For example, we could proceed as follows: If g is equal to the bth nucleobase of the genome among {G anc , G (1) , G (2) } that is the nearest in time from G node . This could lead to a clear improvement when G anc is S exo , because the date t exo of S exo is often far away from dates of The approach presented in this article for reconstructing the virus sequences at the internal nodes (and its possible improvement described above) exploits the knowledge of the transmission tree J (more exactly the state of J along the MCMC algorithm) and the knowledge of a common ancestral sequence S exo . The knowledge of J implies that the topology of the evolutionary tree is known (at each stage of the MCMC algorithm) and, consequently, that we have to solve the so-called small parsimony problem (i.e. finding the most parsimonious labeling of the internal vertices in an evolutionary tree; Jones and Pevzner, 2004, chap. 10). This problem is classically solved by the algorithms of [START_REF] Fitch | Toward defining the course of evolution: minimum change for a specific tree topology[END_REF], [START_REF] Hartigan | Minimum mutation fits to a given tree[END_REF] and [START_REF] Sankoff | Locating the vertices of a steiner tree in an arbitrary metric space[END_REF] that consist of two stages: (1) a bottom-up phase where candidate states for the internal nodes are determined and eventually associated with a score; and (2) a top-down refinement where one makes a choice for the sequence state at the tree root and one decides between the candidate states for the other internal nodes by minimizing a total score. The additional knowledge of a common ancestral sequence S exo allows us to simply perform a top-down phase for reconstructing the internal nodes. Thus, Algorithm 2 that is conditional on J and S exo allows us to divide the number of operations at least by a factor 2 compared with the algorithms mentioned above. In practice, S exo will often have to be preliminary reconstructed, but this will be done once and not at each update in the MCMC algorithm that lead to the modification of the genetic likelihood. It has to be noted that Algorithm 2 is not demonstrated to yield the most parsimonious reconstruction of transmitted sequences, but is expected to yield a relatively parsimonious reconstruction because it is based on successive partial reconstructions (i.e. successive calls of Algorithm 1) and each of these partial reconstructions is the most parsimonious.

The reconstruction of evolutionary trees and ancestral sequences has also been tackled with likelihood-based methods such as the maximum likelihood, empirical Bayes and hierarchical Bayes approaches [START_REF] Hanson-Smith | Robustness of ancestral sequence reconstruction to phylogenetic uncertainty[END_REF][START_REF] Huelsenbeck | Empirical and hierarchical bayesian estimation of ancestral states[END_REF][START_REF] Pirie | Model uncertainty in ancestral area reconstruction: a parsimonious solution?[END_REF]. Despite several similarities between (i) the construction of the likelihood defined in these approaches and (ii) the construction of our likelihood, these likelihoods are different because we make specific assumptions about the dependencies between temporal, spatial and genetic variables. However, from an algorithmic point of view, the numerical approaches developed for the reconstruction of evolutionary trees and ancestral sequences could be adapted to our setting. Thus, we could replace our mixed likelihood/parsimony-based inference algorithm by a fully likelihood-based algorithm. Future studies could be carried out in this direction to compare the performance of these algorithms.

As we mentioned in the introduction, this article gives probabilists and statisticians an invitation for enriching the class of space-time-genetic SEIR models and developing efficient inference algorithms in order to carry out high-impact analyses of infectious diseases due to viruses. This invitation is particularly relevant because, today, (i) our connected world favors in some ways the spread of infectious diseases, and (ii) sequencing genetic material has become so affordable that they are available for many infectious diseases. New statistical methods combining spatial, temporal and genetic data and estimating transmission trees will become very important for future epidemiology. 0.05 0.10 0.15 0.20 0.25 0.30 0.00 0.02 0.04 0.06 0.08 0.10 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q FIGURE 4.

  where γ(•; a, b) is the probability distribution function of the gamma distribution parameterized by (a, b).

FIGURE 3 .

 3 FIGURE 3. Example of transmission tree with five individuals i, j, k, l and m, where j and m are infected by the exogenous source and i, k and l are infected by observed individuals. On horizontal lines, that are time lines, the first tick indicates the infection time, the second the end of exposed stage, and the third the observation time. A: Evolution of sequences are independent on tree segments marked with double arrows; the sequences S in f at infection times appear in the genetic likelihood given in Subsection 3.4 and in its approximation given in Subsection 3.6. B: In the genetic pseudo-likelihood of Subsection 3.5, the probability of observing S in f i , for example, is assessed by computing the probability of evolution between S in f i and S in f k and the probability of evolution between S in f i and S in f j , without taking into account the dependence due to the overlapping of these two evolutions.

TABLE 1 .

 1 Possible events and corresponding transition rates for the discrete-state, continuous-time Markovian SEIR model.

TABLE 2 .

 2 Possible

events and corresponding transition rates for the spatial, discrete-state, continuous-time Markovian SEIR model.

TABLE 3 .

 3 Possible events and corresponding transition rates for the individual-based, spatial, discrete-state, continuous-time Markovian SEIR model.

	Description	Event	Rate
	Infection	S k

TABLE 4 .

 4 Possible events and corresponding transition rates or distributions for the individual-based, spatial, discrete-state, continuous-time, semi-Markov SEIR model.

	Description	Event	Rate	Distribution
	Infection	S k

  δ 2 )

	Soumis au Journal de la Société Française de Statistique
	File: seir.tex, compiled with jsfds, version : 2009/12/09
	date: September 29, 2015

TABLE 5 .

 5 Possible events and corresponding substitution rates for the Markovian evolutionary model. Letters A, C, G and U denotes nucleobases adenine, cytosine, guanine and thymine, respectively. Transversion of type 1 A→U or U→A or C→G or G→C µ 2 Transversion of type 2 A→C or C→A or G→U or U→G µ 3 -The disease reservoir (i.e. exogenous source) is assumed to simply consist of one virus sequence S exo dated at time t exo ∈ R. -We assume that, at any time, there is only one sequence of the virus (i.e. one viral variant) per infected individual. The sequence in individual k at time t, for t such that E k

	Description	Event	Rate
	Transition	A→G or G→A or C→U or U→C µ 1

  Soumis au Journal de la Société Française de Statistique File: seir.tex, compiled with jsfds, version : 2009/12/09 date: September 29, 2015 fraction of infected individuals are observed. Now, consider Figure

  (T in f k ) (the dimension of these unknowns is n × s),[START_REF] Morelli | A Bayesian inference framework to reconstruct transmission trees using epidemiological and genetic data[END_REF] andMollentze et al. (

Soumis au Journal de la Société Française de Statistique File: seir.tex, compiled with jsfds, version : 2009/12/09 date: September 29, 2015 to handle latent vectors S k

  Soumis au Journal de la Société Française de Statistique File: seir.tex, compiled with jsfds, version : 2009/12/09 date: September 29, 2015

  Soumis au Journal de la Société Française de Statistique File: seir.tex, compiled with jsfds, version : 2009/12/09 date: September 29, 2015

  Soumis au Journal de la Société Française de Statistique File: seir.tex, compiled with jsfds, version : 2009/12/09 date: September 29, 2015

	TABLE

  of the set of latent variables lead to inference difficulties. The main difficulty concerns the manner to handle latent genetic sequences (i.e. the transmitted pathogen sequences at infection times). We proposed a new algorithm including a new way to deal with latent genetic sequences.

		-likelihood		Approximate likelihood
	Genetic sampling effort (%)	100	100	50	25
	C 1 (true endogenous transmissions)	0.75 (0.20)	0.82 (0.20) 0.48 (0.27) 0.36 (0.31)
	C 2 (true exogenous transmissions)	0.80 (0.18)	0.80 (0.18) 0.72 (0.21) 0.72 (0.21)
	C 3 (false exogenous transmissions)	0.14 (0.18)	0.11 (0.17) 0.26 (0.26) 0.34 (0.33)

dimension

  Soumis au Journal de la Société Française de Statistique File: seir.tex, compiled with jsfds, version : 2009/12/09 date: September 29, 2015 bifurcation nodes coinciding with infection times of observed individuals.

  Example of simulated epidemic. Black dots represent infected hosts, while black segments represent transmissions. Samples are only taken from a sub-region, colored in gray, and in this area, not all cases are detected or sampled. Infected hosts have a probability of 3/4 of being sampled. Points corresponding to sampled hosts are surrounded by red circles. Soumis au Journal de la Société Française de Statistique File: seir.tex, compiled with jsfds, version : 2009/12/09 date: September 29, 2015

Soumis au Journal de la Société Française de Statistique File: seir.tex, compiled with jsfds, version : 2009/12/09 date: September 29, 2015

Acknowledgements

A large part of the content in this manuscript is a probabilistic and statistical formalization of previous articles, namely Morelli et al. (2012) and Mollentze et al. (2014), and the author wish to thanks all co-authors of these articles. This work was funded by the ModEE grant (SPE division of INRA).

Appendix A: On the distribution of susceptible state durations

In each SEIR model presented in Section 2, the transition rate governing the changes of the susceptible state S → S -1 (resp. S k → S k -1) is a function of the infectious state(s) I (resp. {I j : j = k}). Suppose that these infectious states are constant in time (which is not true in general), then the duration required for a change in S or S k is exponentially distributed. However, I and {I j : j = k} are randomly varying. Therefore, (i) the transition rate governing the changes of the susceptible state S or S k is a random time-varying function and (ii) the duration required for a change in S or S k is not exponential in general.

For instance, consider the stochastic time process {S(t) : t ≥ 0} defined in Subsection 2.1. The times at which S changes (i.e. makes a negative jump of one unit) form a Cox process (or doubly stochastic point process; [START_REF] Daley | An introduction to the theory of point processes, volume i: Elementary theory and methods of probability and its applications[END_REF][START_REF] Grandell | Doubly stochastic Poisson processes[END_REF] where the stochastic intensity function is t → α 0 S(t) + α 1 S(t)I(t)/T(t). Note that this intensity function is zero once S(t) = 0. Closed (and non-exponential) forms of the distribution of inter-point times have only been obtained for specific Cox processes (e.g. [START_REF] Dassios | The distribution of the interval between events of a cox process with shot noise intensity[END_REF]. In the general case, computing the distribution of inter-point times requires an integration over the possible realizations of the intensity function.