Keywords: polygenic mutation, stabilizing selection, mixed mating, inbreeding depression, purging

ABSTRACT We analyze two models of the maintenance of quantitative genetic variance in a mixed mating system of self-fertilization and outcrossing. In both models purely additive genetic variance is maintained by mutation and recombination under stabilizing selection on the phenotype of one or more quantitative characters. The Gaussian Allele Model (GAM) involves a finite number of unlinked loci in an infinitely large population, assuming a normal distribution of allelic effects at each locus within lineages selfed for τ consecutive generations since their last outcross. The Infinitesimal Finite-population Model (IFM) involves an infinite number of loci in a large but finite population, assuming a normal distribution of breeding values in lineages of selfing age τ , with no assumption on the distribution of allelic effects within loci. In both models a stable equilibrium genetic variance exists, the outcrossed equilibrium, nearly equal to that under random mating, for all selfing rates, r, up to critical value, r, the purging threshold, which approximately equals the mean fitness under random mating relative to that under complete selfing. In the GAM a second stable equilibrium, the purged equilibrium, exists for any positive selfing rate, with genetic variance less than or equal to that under pure selfing; as r increases above r the outcrossed equilibrium collapses sharply to the purged equilibrium genetic variance. In the IFM a single stable equilibrium genetic variance exists at each selfing rate; as r increases above r the equilibrium genetic variance drops sharply and then declines gradually to that maintained under complete selfing. We discuss implications of the models for the evolution of selfing rates, and their relevance to the classical view of Stebbins that predominant selfing constitutes an "evolutionary dead end".

Many species of flowering plants, and some hermaphroditic animals, reproduce by a mixture of self-fertilization and outcrossing, often evolving to predominant selfing [START_REF] Stebbins | Self fertilization and population variability in higher plants[END_REF][START_REF] Stebbins | Flowering Plants. Evolution Above the Species Level[END_REF][START_REF] Harder | Ecology and Evolution of Flowers[END_REF][START_REF] Igic | Bias in the studies of outcrossing rate distributions[END_REF][START_REF] Jarne | Animals mix it up too: the distribution of selffertilization among hermaphroditic animals[END_REF]. In such mixed mating systems inbreeding depression for fitness is a critical determinant of mating system evolution [START_REF] Lande | The evolution of self-fertilization and inbreeding depression in plants. I. Genetic models[END_REF][START_REF] Charlesworth | Inbreeding depression and its evolutionary consequences[END_REF][START_REF] Charlesworth | Inbreeding depression, genetic load, and the evolution of outcrossing rates in a multilocus system with no linkage[END_REF][START_REF] Charlesworth | The genetics of inbreeding depression[END_REF][START_REF] Porcher | Loss of gametophytic self-incompatibility with evo-lution of inbreeding depression[END_REF][START_REF] Devaux | Pollination ecology and inbreeding depression control individual flowering phenologies and mixed mating[END_REF]. Spontaneous deleterious mutations, as well as standing genetic variation, contribute to inbreeding depression and display a strongly bimodal distribution of allelic effects on fitness [START_REF] Dobzhansky | Genetics of the Evolutionary Process[END_REF][START_REF] Fudala | Low frequency of mutations with strongly deleterious but nonlethal fitness effects[END_REF][START_REF] Bell | Experimental genomics of fitness in yeast[END_REF]. Lethal and semi-lethal mutations in standing variation are on average nearly recessive, whereas mildly deleterious mutations have slightly recessive to nearly additive fitness effects [START_REF] Simmons | Mutations affecting fitness in Drosophila populations[END_REF][START_REF] Willis | Inbreeding load, average dominance and the mutation rate for mildly deleterious alleles in Mimulus guttatus[END_REF][START_REF] Vassilieva | The fitness effects of spontaneous mutations in Caenorhabditis elegans[END_REF][START_REF] Eyre-Walker | The distribution of fitness effects of new mutations[END_REF].

Stabilizing selection on quantitative characters is thought to be prevalent in natural populations and, although it may be weak or fluctuating on many characters [START_REF] Wright | Evolution and the Genetics of Populations[END_REF][START_REF] Lande | The role of genetic variability in adaptation and population persistence in a changing environment[END_REF][START_REF] Kingsolver | The strength of phenotypic selection in natural populations[END_REF][START_REF] Lande | Expected relative fitness and the adaptive topography of fluctuating selection[END_REF], it may create a substantial component of the total inbreeding depression for fitness. Stabilizing selection on a quantitative character produces allelic effects on fitness that are mildly deleterious and slightly recessive [START_REF] Wright | The analysis of variance and the correlations between relatives with respect to deviations from an optimum[END_REF], in agreement with observations on the mildly deleterious component of inbreeding depression [START_REF] Simmons | Mutations affecting fitness in Drosophila populations[END_REF][START_REF] Willis | Inbreeding load, average dominance and the mutation rate for mildly deleterious alleles in Mimulus guttatus[END_REF]).

Under stabilizing selection on a quantitative character, an allele with an additive effect on the character may be either advantageous or deleterious depending on whether the mean phenotype is above or below the optimum, and alleles at different loci with opposite effects on the character may compensate for each other in their effects on phenotype and fitness, even when the mean phenotype is at the optimum [START_REF] Fisher | The Genetical Theory of Natural Selection[END_REF](Fisher , 1958;;[START_REF] Wright | Evolution in Mendelian populations[END_REF][START_REF] Wright | Evolution and the Genetics of Populations[END_REF]. [START_REF] Charlesworth | Stabilizing selection, purifying selection and mutational bias in finite populations[END_REF] highlighted the difficulty of empirically distinguishing the relative contributions of unconditionally versus conditionally deleterious mutations of small effect.

With respect to primary characters of morphology, physiology and behavior that determine fitness, [START_REF] Wright | Systems of mating. II. The effects of inbreeding on the genetic composition of a population[END_REF][START_REF] Wright | Evolution and the Genetics of Populations[END_REF] showed for purely additive genetic variance in the absence of mutation and selection that inbreeding can increase the genetic variance up to a factor of 2. [START_REF] Lande | The influence of the mating system on the maintenance of genetic variability in polygenic characters[END_REF] modelled the maintenance of quantitative genetic variance by mutation under stabilizing selection for regular systems of non-random mating, including inbreeding with no variance in inbreeding coefficient among individuals. In sharp contrast with Wright's classical results, [START_REF] Lande | The influence of the mating system on the maintenance of genetic variability in polygenic characters[END_REF] found that a regular system of non-random mating has no impact on the equilibrium genetic variance maintained by mutation and stabilizing selection. [START_REF] Charlesworth | Quantitative genetics in plants: the effect of the breeding system on genetic variability[END_REF] modelled the maintenance of quantitative variation by mutation under random mating versus complete selfing for different models of selection, concluding that complete selfing substantially reduces the genetic variance compared to random mating, but they did not model mixed mating.

Mixed mating systems, such as combined selfing and outcrossing, introduce the serious complication of zygotic disequilibrium (non-random associations of diploid genotypes among loci) and variance of inbreeding coefficient among individuals [START_REF] Haldane | The association of characters as a result of inbreeding and linkage[END_REF][START_REF] Crow | An Introduction to Population Genetics Theory[END_REF]. These complications are successfully encompassed only by models of unconditionally deleterious mutations (Kondrashov 1985, Charlesworth and[START_REF] Charlesworth | Some evolutionary consequences of deleterious mutations[END_REF]. A theory of the maintenance of quantitative genetic variance does not currently exist for such mixed mating systems. To account for zygotic disequilibrium in quantitative characters under mixed mating, we employ the selfing age structure of the population, introduced by [START_REF] Kelly | Mutation-selection balance in mixed mating populations[END_REF] into the [START_REF] Kondrashov | Deleterious mutations as an evolutionary factor. II. Facultative apomixis and selfing[END_REF] model.

We analyze two models of the maintenance of quantitative genetic variance in a mixed mating system of self-fertilization and outcrossing. In both models purely additive genetic variance is maintained by mutation and recombination among unlinked loci under stabilizing selection on the phenotype of one or more quantitative characters. The Gaussian Allele Model (GAM) involves a finite number of unlinked loci in an infinitely large population, assuming a normal distribution of allelic effects at each locus within lineages selfed for τ consecutive generations since their last outcross. The Infinitesimal Finite-population Model (IFM) involves an infinite number of loci in a very large but finite population, assuming a normal distribution of breeding values in lineages of selfing age τ , with no assumption on the distribution of allelic effects within loci. Aspects of the results common to both models are considered to be robust, and have fundamental implications for the evolution of plant mating systems.

Basic assumptions of the models

The diploid population is partially self-fertilizing such that each zygote has a probability r of being produced by self-fertilization and probability 1 -r of being produced by outcrossing to an unrelated individual. There is no genetic variation in the selfing rate. Genetic variance in quantitative traits under stabilizing selection is assumed to be purely additive. The population is measured in each generation before selection on a (vector of) quantitative character(s), z. Under random mating the total additive genetic variance (or variance-covariance matrix) can be partitioned into two additive components G = V + C. Diagonal elements of V give the genic variances of each character (twice the variance of allelic effects on each character summed over all loci), and off-diagonal elements of V give the genetic covariances between characters due to pleiotropy (twice the covariances of allelic effects between pairs of characters summed over all loci with pleiotropic effects on the characters). C is the variance-covariance matrix of twice the total covariance of allelic effects among loci within gametes due to linkage disequilibrium (nonrandom association of alleles between loci within gametes). In the absence of selection, inbreeding reduces heterozygosity, proportionally reducing additive genetic (co)variance within families and increasing additive genetic (co)variance among families [START_REF] Wright | Systems of mating. II. The effects of inbreeding on the genetic composition of a population[END_REF][START_REF] Wright | Evolution and the Genetics of Populations[END_REF][START_REF] Crow | An Introduction to Population Genetics Theory[END_REF].

Individual environmental effects on the phenotype are assumed to be independent among individuals, normally distributed with mean 0 and variance E, and additive and independent of the selfing age or breeding value (total additive genetic effect summed over all loci in an individual). In any generation before phenotypic selection a cohort of selfing age τ has genetic and phenotypic variances respectively of

G τ = (1 + f τ )(V τ + C τ ) (1a)
P τ = G τ + E (1b)
Here f τ is [START_REF] Wright | Systems of mating. II. The effects of inbreeding on the genetic composition of a population[END_REF][START_REF] Wright | Evolution and the Genetics of Populations[END_REF] biometrical correlation of additive effects of alleles at the same locus within individuals, rather than Malécot's probability of identity by descent. These two measures of f may be rather different since small mutational changes in additive effects of alleles cause only a small decrease in the biometrical correlation between alleles at the same locus, whereas mutation completely eliminates allelic identity. Assumptions specific to each model guarantee that for all individuals within a cohort of selfing age τ a single value of f τ applies to all loci affecting a given character (GAM) or to all characters (IFM). The covariance of additive effects of alleles from different gametes equals f τ multiplied by the variance of allelic effects within gametes (or the covariance of allelic effects at different loci within gametes).

The genetic and phenotypic variances before selection in the population as a whole are denoted without subscripts as

G = ∞ τ =0 p τ G τ (1c)
P = G + E (1d)
where p τ is the frequency of selfing age τ in the population before selection.

Multivariate stabilizing selection on the individual phenotype, z, is described by a Gaussian function of the individual deviation from an optimum phenotype, θ,

W (z) = exp -1 2 (z -θ) T Ω -1 (z -θ) (2)
where Ω is a symmetric matrix describing stabilizing and correlational selection [START_REF] Lande | The measurement of selection on correlated characters[END_REF], with where superscripts T and -1 respectively denote vector or matrix transpose and matrix inverse. Using the normal approximation for phenotypes and breeding values, the mean fitness of a cohort of selfing age τ is then

wτ = φ τ (z)W (z)dz = |Ωγ τ | exp -1 2 (z -θ) T γ τ (z -θ)
where φ τ (z) is the phenotypic distribution in the cohort of selfing age τ . Vertical bars | | denote the determinant of a matrix and γ τ = (Ω + P τ ) -1 . In subsequent formulas diagonal elements of the Ω matrix are denoted as ω 2 , and for independently selected characters the off-diagonal elements are 0. The general environment is assumed to be constant such that the mean phenotype always remains at the optimum, z = θ, so the mean fitness of a selfing age cohort simplifies to

wτ = |Ωγ τ | (3a)
and the population mean fitness is

w = ∞ τ =0 p τ wτ (3b)
For independently selected characters, the mean fitness within each selfing age class is the product of the mean fitnesses of the characters, but this is not true for the population as a whole.

Gaussian Allele Model (GAM)

Mutation-selection balance for one character

We employ the mutation model of [START_REF] Kimura | A stochastic model concerning the maintenance of genetic variability in quantitative characters[END_REF] and [START_REF] Lande | The maintenance of genetic variability by mutation in a polygenic character with linked loci[END_REF][START_REF] Lande | The influence of the mating system on the maintenance of genetic variability in polygenic characters[END_REF] with n loci having completely additive effects on a particular character. Each allele at a given locus mutates at a constant rate with the same distribution of changes in additive effect on the character, with no directional bias. This produces a constant mutational variance for the character, σ 2 m , without changing the mean phenotype. Empirical estimates of the mutational variance, scaled by the environmental variance (or as here using E = 1) are typically about 10 -3 to 10 -4 [START_REF] Lande | The maintenance of genetic variability by mutation in a polygenic character with linked loci[END_REF][START_REF] Charlesworth | Quantitative genetics in plants: the effect of the breeding system on genetic variability[END_REF][START_REF] Lynch | The rate of polygenic mutation[END_REF][START_REF] Houle | Comparing mutational variabilities[END_REF].

Assuming weak stabilizing selection and a high mutation rate per locus, the distribution of allelic effects at each locus is approximately Gaussian [START_REF] Kimura | A stochastic model concerning the maintenance of genetic variability in quantitative characters[END_REF][START_REF] Bürger | The Mathematical Theory of Selection, Recombination and Mutation[END_REF]. This is consistent with empirical observations of high genomic mutation rates per character, on the order of 10 -1 to 10 -2 observed for quantitative traits in maize [START_REF] Russell | Mutations affecting quantitative characters in long-time inbred lines of maize[END_REF], and the assumption that n is much less than the total number of genes in the genome, on the order of n = 10 to 100 for the effective number of loci [START_REF] Lande | The maintenance of genetic variability by mutation in a polygenic character with linked loci[END_REF][START_REF] Lande | The influence of the mating system on the maintenance of genetic variability in polygenic characters[END_REF].

For simplicity, we assume n unlinked loci with equal mutational variance, so that a single inbreeding coefficient applies to all loci affecting a given character within a given age class. Multiple characters are assumed to be genetically and phenotypically independent and subject to independent stabilizing selection. For multiple characters that differ in their parameters (number of loci, mutational variance, and stabilizing selection) a different set of recursions across the selfing age classes is required.

Gamete production from the selfing age classes

Summing the additive effects of alleles at exchangeable loci in [START_REF] Lande | The maintenance of genetic variability by mutation in a polygenic character with linked loci[END_REF][START_REF] Lande | The influence of the mating system on the maintenance of genetic variability in polygenic characters[END_REF], the genetic covariance c τ and genic variance v τ of gametes produced by individuals in selfing age class τ ≥ 1 are

c τ = 1 2 1 + f τ 2 C τ -1 - 1 n γ τ 2 G 2 τ (4a) v τ = 1 2 V τ - γ τ 2n G 2 τ + σ 2 m (4b)
For the outcrossed class τ = 0, a straightforward way to derive the components of genetic variance and gametic output is through a weighted average of gametes produced by all selfing age classes. However, this approach mixes distributions with possibly Outcrossing by random mating among selfing age classes implies that for pairs of parents with selfing ages i and j the frequency of the ij subclass of the outcrossed class 0, denoted as p 0ij , is simply the product of the frequencies of parental age classes at the adult stage, after selection in the previous generation denoted by subscript (t -1),

p 0ij = wi(t-1) w(t-1) p i(t-1) wj(t-1) w(t-1) p j(t-1) (5)
The mean fitness of the outcrossed class as a whole is then

w0 = ∞ i=0 ∞ j=0 p 0ij w0ij (6)
The mean fitness of each outcrossed subclass, w0ij , depends only on its total genetic variance, but selection acts differently on unequal gametic contributions to the genic variance and covariance by parents of different selfing age. Gametes produced by the outcross class, averaged over all subclasses, have genic variance and covariance

c 0 = 1 2 w0 ∞ i=0 ∞ j=0 p 0ij w0ij c i(t-1) + c j(t-1) 2 (7a) -γ 0ij 1 - 1 n [c i(t-1) + v i(t-1) ] 2 + [c j(t-1) + v j(t-1) ] 2 v 0 = 1 2 w0 ∞ i=0 ∞ j=0 p 0ij w0ij v i(t-1) + v j(t-1) - γ 0ij n [c i(t-1) + v i(t-1) ] 2 + [c j(t-1) + v j(t-1) ] 2 + σ 2 m 2 (7b)
where subscript (t -1) denotes the grandparental generation. The relative frequencies and mean fitness of class 0ij are obtained using the components of genetic variance for the outcrossed progeny of a mating between individuals of selfing age classes i and j, with a prime denoting the next generation,

C 0ij = c i + c j (8a) V 0ij = v i + v j (8b) f 0ij = 0. ( 8c
)
These yield the total genetic and phenotypic variance (using eqs. 1), and hence the mean fitness w0ij (eqs. 3) of the subclass.

Recursions for components of genetic variance

All selfing age classes after the first obey the recursions, for τ ≥ 1,

C τ +1 = 2c τ = 1 + f τ 2 C τ -1 - 1 n γ τ 2 G 2 τ (9a) V τ +1 = 2v τ = V τ - γ τ 2n G 2 τ + σ 2 m (9b) f τ +1 = 1 2v τ 1 + f τ 2 V τ - γ τ 2n G 2 τ (9c)
In the absence of selection and mutation eq. ( 9c) reduces to the recursion of [START_REF] Wright | Systems of mating. II. The effects of inbreeding on the genetic composition of a population[END_REF][START_REF] Wright | Evolution and the Genetics of Populations[END_REF] for the inbreeding coefficient under continued selfing,

f τ +1 = (1 + f τ )/2.
The genetic variance components of the outcrossed class are twice the weighted average of gametic outputs from all selfing age classes:

C 0 = 2 ∞ τ =0 p τ wτ w c τ (10a) V 0 = 2 ∞ τ =0 p τ wτ w v τ (10b) f 0 = 0 (10c)
Finally, for the first selfing age class τ = 1,

C 1 = 2c 0 (11a) V 1 = 2v 0 (11b) f 1 = 1 2v 0 ∞ i=0 ∞ j=0 p 0ij w0ij w0 v i(t-1) + v j(t-1) 2 - γ 0ij n [c i(t-1) + v i(t-1) ] 2 + [c j(t-1) + v j(t-1) ] 2 (11c)

Age distribution of selfing lineages

After stabilizing selection, mating and reproduction, the distribution of selfing ages in the population is

p 0 = 1 -r (12a) p τ +1 = r wτ w p τ (12b)
Eqs. ( 1)-( 12) constitute the complete recursion system for the evolution of quantitative genetic variance. For numerical computation it is necessary to truncate the age distribution of selfing lineages at some upper limit. This approach is often used in analyzing the demography of age-structured populations [START_REF] Caswell | Matrix Population Models[END_REF]. Recursion formulas for the frequency and genetic composition of the terminal selfing age class are

given in Appendix 1. The number of classes retained for numerical analysis should be sufficiently large that substantially increasing it does not appreciably affect the results.

Numerical analysis of the GAM for multiple genetic independent characters under independent stabilizing selection showed that the results had a low sensitivity to substantial differences in stabilizing selection or mutational variance among characters (unpublished results).

Infinitesimal Finite-Population Model (IFM)

To relax the critical but controversial assumption in the GAM of a normal distribution of allelic effects within loci [START_REF] Turelli | Heritable genetic variation via mutation-selection balance: Lerch's zeta meets the abdominal bristle[END_REF][START_REF] Turelli | Genetic and statistical analyses of strong selection on polygenic traits: What, me normal?[END_REF], and to facilitate analysis of correlated characters, we extend [START_REF] Fisher | The correlation between relatives on the supposition of Mendelian inheritance[END_REF] An unrealistic feature of the classical infinitesimal model of [START_REF] Fisher | The correlation between relatives on the supposition of Mendelian inheritance[END_REF] and [START_REF] Bulmer | The effect of selection on genetic variability[END_REF] for an infinite population is that the equilibrium inbreeding depression due to stabilizing selection increases with the selfing rate of a population, until at sufficiently high selfing rates stabilizing selection finally creates enough negative linkage disequilibrium to purge the genetic variance (unpublished results). This unrealistic feature of the classical infinitesimal model for an infinite population can be understood from the classical result of [START_REF] Wright | Systems of mating. II. The effects of inbreeding on the genetic composition of a population[END_REF][START_REF] Wright | Evolution and the Genetics of Populations[END_REF] for an infinite population with no selection or mutation, in which the equilibrium genetic variance increases as a linear function of the population inbreeding coefficient.

To obtain a more realistic infinitesimal model for partial selfing, we introduce the Infinitesimal Finite-population Model (IFM), in which the genic variance, V , is maintained by a balance between mutation and random genetic drift. The IFM can then incorporate the well-known influence of inbreeding in changing the effective population size and hence the genic variance maintained by mutation.

The IFM is derived from the GAM by letting the number of loci approach infinity, n → ∞, and, simultaneously letting the effective population size under random mating become very large and the mutational variance become very small, N e (0) → ∞ and σ 2 m → 0, such that N e (0)σ 2 m remains constant. Under random mating the IFM has the same dynamics in response to selection as in the classical infinitesimal model of [START_REF] Fisher | The correlation between relatives on the supposition of Mendelian inheritance[END_REF] and [START_REF] Bulmer | The effect of selection on genetic variability[END_REF], but more generally the IFM also allows the genic variance to adjust to changes in the selfing rate as follows. With these assumptions the recursions for the inbreeding coefficient as a function of selfing age, eqs (9c), (10c) and (11c) reduce to Wright's classical formula for continued selfing

f τ +1 = (1 + f τ )/2 (13a)
The genic variance in the total population (or in any selfing age class) obeys the recursion

V (r) = 1 - 1 2N e (r) V (r) + σ 2 m
where N e (r) is the effective population size at selfing rate r. It is well-known that a completely selfing population has an effective size half that under random mating [START_REF] Charlesworth | Quantitative genetics in plants: the effect of the breeding system on genetic variability[END_REF]. More generally, using methods of [START_REF] Wright | Evolution in Mendelian populations[END_REF][START_REF] Wright | Evolution and the Genetics of Populations[END_REF], the effective size of a population with inbreeding coefficient f is N e (r) = N e (0)/(1 + f ). Substituting this into the recursion for the genic variance and using [START_REF] Wright | Systems of mating. II. The effects of inbreeding on the genetic composition of a population[END_REF][START_REF] Wright | Evolution and the Genetics of Populations[END_REF] formula for the equilibrium inbreeding coefficient in a partial selfing population with no selection or mutation, f = r 2-r , gives, asymptotically, the equilibrium genic variance maintained by mutation-drift balance in a large partially selfing population,

V (r) = (1 -r/2)V (0) (13b)
where V (0) = 2N e (0)σ 2 m in agreement with previous results on mutation-selection balance under random mating [START_REF] Clayton | Mutation and quantitative variation[END_REF][START_REF] Lande | Genetic variation and phenotypic evolution during allopatric speciation[END_REF]. For any given selfing rate, this equilibrium genic variance replaces the recursions (9b), ( 10b) and (11b). It is half as large for complete selfing as under random mating. Formula (13b), with eqs. ( 1a) and (1c) and Wright's relation f = r 2-r , implies that at equilibrium in the absence of selection (for which C = 0) the genetic variance maintained in the population is actually independent of the selfing rate, G(r) = V (0). This occurs because the reduction of genic variance with larger selfing rate (and smaller effective population size) is exactly compensated by the increased inbreeding coefficient in the population. As we show below, this allows stabilizing selection to purge quantitative genetic variance from the population at high selfing rates, and produces a realistic equilibrium inbreeding depression which always decreases with increased selfing rate.

A great advantage of the IFM is that it can readily model the maintenance of genetic variability in correlated characters under multivariate selection, with genetic correlations between characters due to a combination of pleiotropic mutation and correlational selection. A key ingredient of this generality is that in the IFM Wright's recursion for the inbreeding coefficient under continued selfing (eq. 13a) applies to all loci in the genome regardless of their pleiotropic effects on different characters.

The IFM was used to investigate how pleiotropic mutation and correlational selection changed the results. With many loci of low mutability, in the limit as n → ∞ with very large effective population size, such that V (0) remains constant, the GAM can be converted to a general multivariate IFM by interpreting various quantities as matrices rather than scalars, e.g. in eqs. (4) rewriting γ τ G 2 τ as G τ γG τ [START_REF] Lande | Genetic variation and phenotypic evolution during allopatric speciation[END_REF][START_REF] Lande | The measurement of selection on correlated characters[END_REF]). In the IFM, mutational and genic variance-covariance matrices can be simply transformed by rotation of axes to produce either selectively or genically independent characters (using eigenvectors of the correlational selection matrix γ, or the genic variance-covariance matrix V (0). Our numerical analysis of the IFM therefore focused on mutationally and genically independent characters under correlational selection.

In both models deviation from the assumption of a Gaussian distribution for either the allelic effects (GAM) or breeding values (IFM) within selfing age chorts, occurs due to mixing of genetic contributions among all selfing ages upon outcrossing. We considered the models to have good accuracy when the equilibrium kurtosis in the population remained close to that for a normal distribution (κ = 3). For comparison to the numerical results we derived the equilibrium kurtosis of breeding values in [START_REF] Wright | Systems of mating. II. The effects of inbreeding on the genetic composition of a population[END_REF][START_REF] Wright | Evolution and the Genetics of Populations[END_REF]

Analytical and Numerical Results

Purging genetic variance by stabilizing selection under continued selfing

In both the GAM and the IFM under continued selfing, stabilizing selection purges the genetic variance, but the details of how this happens and the extent of the purging differ in the two models, as shown below.

Continued selfing

In the GAM, assuming weak selection (ω 2 P τ so that γ τ ≈ 1/ω 2 ) and small mutational variance (σ 2 m E) eq. ( 9c) can be expanded as a Taylor series to first order in small terms and solved for a slowly changing quasi-equilibrium, which gives

f τ ≈ 1 -2σ 2 m /V τ .
This can be used with eqs. (9a), ( 9b) and ( 1) to find a first-order approximation for the asymptotic dynamics of the genetic variance and its components for large selfing age, G τ +1 -G τ ≈ -G 2 τ /ω 2 + 4σ 2 m . Thus under continued selfing in the GAM the genetic variance approaches an equilibrium

G ∞ ≈ 2 σ 2 m ω 2 (14a)
This is smaller by a factor of 2/n than the weak selection approximation for the equilibrium genetic variance under random mating G 0 ≈ 2nσ 2 m ω 2 [START_REF] Kimura | A stochastic model concerning the maintenance of genetic variability in quantitative characters[END_REF][START_REF] Lande | The maintenance of genetic variability by mutation in a polygenic character with linked loci[END_REF][START_REF] Lande | The influence of the mating system on the maintenance of genetic variability in polygenic characters[END_REF].

For the GAM the above results with eq. (9a) show that with continued selfing, although G τ approaches a constant, both V τ and -C τ continue to increase indefinitely at the asymptotic rate

lim τ →∞ (V τ +1 -V τ ) = -lim τ →∞ (C τ +1 -C τ ) = (1 - 2 n )σ 2 m (14b)
In view of the dynamics of f τ above eq. ( 14a), this confirms that with increasing selfing age the inbreeding coefficient approaches 1.

For the IFM, a similar analysis produces the asymptotic recursion for the genetic variance under continued selfing, G τ +1 -G τ ≈ -G 2 τ /ω 2 , the only solution of which is

G ∞ = 0. ( 15
)
Under continued selfing, purging of genetic variance in the IFM occurs solely by the buildup of negative linkage from stabilizing selection; as f τ → 1, increasing homozygosity reduces the effective recombination in proportion to 1 -f τ [START_REF] Lande | The influence of the mating system on the maintenance of genetic variability in polygenic characters[END_REF]) so that C τ → -V . By comparison, under random mating in the IFM, assuming weak stabilizing selection (γV 1), the equilibrium genetic variance is approximately

G 0 ≈ (1 -γV )V .
These analytical results concerning the genetic variance and its components as functions of selfing age were confirmed numerically. Fig. 1A shows for the GAM the indefinite increase of the genic variance and the negative linkage disequilibrium variance with increasing selfing age, due to the buildup of associations between compensatory mutations (eq. 14b) while the genetic variance approaches a constant (eq. 14a). For populations with a high selfing rate the accumulation of a large negative linkage disequilibrium in the GAM leads to recombination and segregation in the second generation after outcrossing (τ = 1) producing a high genetic variance and low mean fitness (Fig. 1B,D). Similar effects occur in the IFM, but with restricted magnitude (Fig. 2B,D).

[Fig. 1 and Fig. 2 about here]

Mixed mating and the distribution of selfing ages

With partial selfing, the distribution of selfing ages in the population plays an crucial role in the dynamics of genetic variance. At a stable equilibrium the distribution of selfing ages is p τ = (1 -r)r τ l τ where (1 -r)r τ is the probability of selfing τ generations in a row, and l τ = w-τ τ -1 i=0 wi for τ ≥ 1 is the relative survival of selfed lineages to selfing age τ , with l 0 = 1.

Despite continued selfing and stabilizing selection eventually purging the genetic variance in both models, at any given selfing rate in the population, the distribution of selfing ages determines the extent of purging of the genetic variance in the population as a whole. This occurs because shifts in the selfing age distribution change the balance of contributions of young and old selfing ages to the genetic variance at outcrossing (eq. 9) which is then transmitted through the selfing ages (eqs. 8, 10). A low or moderate r shifts the selfing age distribution to the left, toward younger ages; a high r shifts the selfing age distribution right, toward older ages. A crucial property of the relative survival function of selfed lineages is that the mean fitness at each selfing age depends on selection on the whole organism rather than just a single character (eq. 11); this explains why selection on multiple characters affects the purging of genetic variance in each character, as illustrated in the numerical results.

Numerical examples of the mean fitness as a function of selfing age, and the selfing age distribution, are illustrated in Fig. 1C-F and Fig. 2C-F for the GAM and the IFM respectively. At high selfing rate in the GAM the long-selfed lineages become reproductively isolated from the outcrossed lineages, as shown by the very low fitness of intermediate selfing ages and the nearly disjunct bimodal distribution of selfing ages (Fig. 1D,F). Similar effects occur in the IFM at high selfing rate, but to a lesser extent (Fig. 2D,F).

Equilibrium genetic variance in the population as a function of r

For intermediate selfing rates the complexity of the models precludes analytical solution. Recursions for the GAM and IFM were iterated numerically to characterize the equilibrium genetic variance in the population as a function of population selfing rate for a wide range of parameter values. Below we describe the salient similarities and differences between the GAM and IFM, emphasizing the robust results common to both models. The accuracy of the assumption of a normal distribution of breeding values within selfing age cohorts in both models was assessed using the kurtosis of breeding values for the population as a whole, κ. To elucidate patterns in the equilibrium genetic variance, we computed for the population as a whole the mean fitness (eq. 3b), the inbreeding depression, δ (loss of mean fitness upon selfing vs. outcrossing), and the genetic variance after selection within a generation, G * (Appendix 2).

In the GAM, two types of stable equilibrium exist for the total genetic variance in the population. Outcrossed equilibrium: Over a wide range of selfing rates from 0 up to a critical selfing rate, r, termed the purging threshold, the genic variance, V , and linkage disequilibrium, C, evolve to nearly compensate for non-random mating, maintaining nearly the same total genetic variance, G, as under random mating. At selfing rates above r the stable outcrossed equilibrium collapses to the purged equilibrium. Purged equilibrium: Selfing rates even slightly above the purging threshold cause a collapse of the equilibrium genetic variance after selection within a generation, G * , to values equal to or less than under pure selfing (eq. 14a). A stable purged equilibrium exists for all selfing rates. The purged equilibrium exists and is stable at any selfing rate (r > 0), but its stability becomes weaker and its domain of attraction smaller for lower selfing rates.

At selfing rates below the purging threshold, r < r, the initial condition of genetic monomorphism always leads to the outcrossed equilibrium; the purged equilibrium is attained from initial conditions with large negative linkage disequilibrium and small total genetic variance. At low selfing rates the rate of convergence to the purged equilibrium may be slow and its domain of attraction small.

[Fig. 3 and Fig. 4 about here] Fig. 3 illustrates for the GAM that at selfing rates below the purging threshold, under weak stabilizing selection the equilibrium genetic variance after selection within a generation, G * , is only slightly smaller than that before selection, G. For selfing rates above the purging threshold, the equilibrium genetic variance before selection, G, blows up to infinity. This occurs because of the unlimited negative linkage disequilibrium (compensatory mutations) and genic variance built up under continued selfing, which segregates in the second generation after outcrossing and recombination. However, the F 2 recombinants among long-selfed lineages are strongly selected against because of their large phenotypic deviations from the optimum. For this reason, at selfing rates above the purging threshold, the equilibrium genetic variance after selection within a generation, G * , collapses to G * purged , the purged equilibrium after selection. For the same reason, at selfing rates above the purging threshold the equilibrium kurtosis of breeding value in the population, κ, blows up, but the kurtosis after selection κ * nearly equals that for Wright's neutral model. Thus, under the basic assumptions of the GAM, we consider the numerical results to be reasonably accurate at all selfing rates. Fig. 4 shows that for the IFM the equilibrium genetic variance in the population, G, also nearly equals that under random mating for selfing rates up to a purging threshold. Just above the purging threshold the equilibrium genetic variance dips sharply and then declines gradually with increasing r. The IFM displays only a single stable equilibrium genetic variance at every selfing rate. The kurtosis of breeding value in the population, both before and after selection, becomes large at selfing rates substantially above the purging threshold. Because excess kurtosis increases the strength of selection on the variance [START_REF] Turelli | Adaptive landscapes, genetic distance and the evolution of quantitative characters[END_REF], the IFM overestimates the genetic variance maintained at selfing rates above the purging threshold. Qualitative differences between the models at high selfing rates arise from constancy of the genic variance at any given selfing rate in the IFM, which limits the buildup of negative linkage disequilibrium under continued selfing (compare Fig. 1A to Fig. 2A).

[Fig. 5 and Fig. 6 about here] Fig. 5 and Fig. 6 depict, respectively for the GAM and IFM, how stabilizing selection on multiple independent characters acts synergistically to reduce the purging threshold, and, especially in the GAM, to sharpen the threshold in comparison to that for a single character under the same intensity of stabilizing selection per character.

These figures also display the equilibrium inbreeding depression, δ, and the mean fitness in the population, as functions of the selfing rate. For both models the mean fitness of the population at high selfing rates exceeds that at low selfing rates, consistent with the purging of genetic variance at high selfing rates. Fig. 7 shows analogous results for the IFM with multiple characters with no pleiotropy but under correlational selection on independent pairs of characters. The genetic covariance between pairs of mutually selected characters, B, is then caused solely by linkage disequilibrium, and at high selfing rates their equilibrium genetic correlation, B/G, makes the multivariate distribution of breeding values conform closely to the shape of the fitness surface.

[Fig. 7 about here]

General approximation for the purging threshold

For selfing rates below the purging threshold, in both the GAM and IFM the equilibrium mean fitness in the population at the outcrossed equilibrium remains nearly the same as under random mating, Wout . More remarkably, in both models at selfing rates above the purging threshold the mean fitness nearly equals the product of the selfing rate, r and the equilibrium mean fitness under completely selfing, Wself . This can be seen most explicitly from the dotted lines for the purged equilibrium in Fig. 5C, and by extrapolation of the corresponding lines in Fig. 6C and Fig. 7D from their values at complete selfing back to the origin. The simplicity of these results indicate that in both models the purging threshold, r can be accurately located by the intersection of these two lines, Wout ≈ r Wself . Remarkably, the purging threshold can thus be accurately approximated as the ratio of mean fitnesses at equilibrium under random mating versus pure selfing,

r ≈ Wout Wself . (16)
With many characters, the numerical analysis indicates a sharp purging threshold, and the analytical approximation is fairly accurate. A small inaccuracy arises, most noticeably in the GAM for m = 5 or 10, because Wout as a function of r dips slightly near the intersection with r Wself .

purged equilibrium becomes reproductively isolated from its outcrossing relatives. This mechanism for speciation by predominant selfing accords with the finding of [START_REF] Goldberg | Tempo and mode in plant breeding system evolution[END_REF] that phylogenetic transitions to predominant selfing often coincide with species originations.

In the IFM only a single stable equilibrium for the genetic variance exists at any selfing rate. This qualitative difference between the two models occurs because in the IFM the genic variance at any given selfing rate remains constant, which limits the buildup of linkage disequilibrium variance under continued selfing; this also limits the magnitude of F 2 breakdown following outcrossing of a long-selfed lineage.

Another feature common to both the GAM and IFM is that the purging threshold approximately equals the ratio of the equilibrium mean fitness under random mating relative to that under pure selfing. In a constant environment stabilizing selection of multiple quantitative characters toward a joint optimum phenotype always produces a higher mean fitness for the population with the lowest genetic variance, which here is that under pure selfing. A further advantage of selfing is that it also increases mean fitness by building adaptive genetic correlations between characters under correlational selection as shown in Fig. 7 (as in Lande 1984 for a symmetric GAM under a regular system of inbreeding). The complexity of the phenotype, and the overall strength of stabilizing selection on it, are thus of paramount importance in determining the purging threshold. For a single character under moderately strong stabilizing selection the purging threshold occurs at a very high selfing rate (Fig. 2A), but for a complex phenotype with many characters under independent stabilizing selection, the purging threshold is substantially reduced (Fig. 7A).

These results are consistent with limited available data on genetic variance in quantitative characters and inbreeding depression as a function of r [START_REF] Charlesworth | Quantitative genetics in plants: the effect of the breeding system on genetic variability[END_REF]. They also have important implications for evolution of the selfing rate, discussed below. In this regard it is important to note that within any very large population different sets of quantitative characters will conform more closely to the assumptions of either the GAM or IFM. Even if only a single character conforms to the GAM, a stable purged equilibrium will exist at all positive selfing rates (Fig. 4A).

However, finite population must eventually limit the buildup of linkage disequilibrium variance; the IFM is therefore likely to be more realistic for most quantitative characters, especially in populations with effective sizes that are not very large.

Genetic variance and inbreeding depression

Inbreeding depression is generally thought to be the main genetic factor counteracting Fisher's automatic advantage of selfing, and preventing the evolution of complete selfing in self-compatible species [START_REF] Fisher | Average excess and average effect of a gene substitution[END_REF][START_REF] Lande | The evolution of self-fertilization and inbreeding depression in plants. I. Genetic models[END_REF][START_REF] Porcher | Loss of gametophytic self-incompatibility with evo-lution of inbreeding depression[END_REF][START_REF] Devaux | Pollination ecology and inbreeding depression control individual flowering phenologies and mixed mating[END_REF]. Previous theory indicates that inbreeding can cause rapid purging of inbreeding depression due to nearly recessive lethals, and that purging of inbreeding depression due to slightly dominant nearly additive mildly deleterious mutations is more difficult, and occurs to a lesser extent [START_REF] Lande | The evolution of self-fertilization and inbreeding depression in plants. I. Genetic models[END_REF][START_REF] Charlesworth | Inbreeding depression, genetic load, and the evolution of outcrossing rates in a multilocus system with no linkage[END_REF][START_REF] Lande | High inbreeding depression, selective interference among loci, and the threshold selfing rate for purging recessive lethal mutations[END_REF][START_REF] Charlesworth | The genetics of inbreeding depression[END_REF][START_REF] Porcher | Evaluating a simple approximation to modeling the joint evolution of self-fertilization and inbreeding depression[END_REF]. Consistent with this is the observation that highly selfing species show substantially reduced inbreeding depression [START_REF] Husband | Evolution of the magnitude and timing of inbreeding depression in plants[END_REF]. [START_REF] Winn | Analysis of inbreeding depression in mixed-mating plants provides evidence for selective interference and stable mixed mating[END_REF] found that species with intermediate selfing rates maintain a substantial inbreeding depression. As they suggest, this may in part be due to selective interference with the purging of lethals that occurs at a high total inbreeding depression [START_REF] Lande | High inbreeding depression, selective interference among loci, and the threshold selfing rate for purging recessive lethal mutations[END_REF], especially if the average inbreeding depression estimated for species with intermediate selfing rates, δ = 0.58, substantially underestimates actual inbreeding depression, due to common omission of the final growth-to-flowering component of the life cycle. [START_REF] Porcher | Evaluating a simple approximation to modeling the joint evolution of self-fertilization and inbreeding depression[END_REF] showed that a larger constant background inbreeding depression produces elevated selective interference among lethals.

The constancy of quantitative genetic variance and inbreeding depression maintained across a wide range of selfing rates below the purging threshold in the present theory is qualitatively consistent with these observations. [START_REF] Charlesworth | Quantitative genetics in plants: the effect of the breeding system on genetic variability[END_REF] reviewed limited available evidence indicating that predominantly selfing populations maintain lower quantitative genetic variance than predominant outcrossers on average. This accords with the reduced genetic variance after selection maintained at selfing rates above the purging threshold in the long-selfed core of the population, and in the population as a whole after selection. In the GAM, at the purged equilibrium after F 2 breakdown to near zero fitness implies that the genetic variance after selection nearly equals the weighted average of that in the long-selfed core and their F 1 outcrossed progeny (Fig. 7A).

In this context it is important to distinguish the purging threshold for recessive lethal mutations [START_REF] Lande | High inbreeding depression, selective interference among loci, and the threshold selfing rate for purging recessive lethal mutations[END_REF]) from that for quantitative genetic variance described here. Although both thresholds involve the reduction of an equilibrium genetic variance at higher selfing rates, they arise under different conditions by different mechanisms, and occur at different selfing rates. Selective interference among lethal mutations depends on a high total inbreeding depression to produce a sharp purging threshold for lethals. In contrast, the purging threshold for quantitative genetic variance exists even under weak stabilizing selection on a single character with consequently small inbreeding depression. The mechanism creating the purging threshold for lethals is that at selfing rates below the threshold nearly all products of selfing are lethal [START_REF] Ganders | Heterozygosity for recessive lethals in homosporous fern populations: Thelypteris palustris and Onoclea sensibilis[END_REF][START_REF] Lande | High inbreeding depression, selective interference among loci, and the threshold selfing rate for purging recessive lethal mutations[END_REF]. In contrast, at selfing rates above the purging threshold for quantitative variance, the only stable state is the purged equilibrium, characterized by F 1 hybrid vigor in crosses between highly inbred lineages at the population's core, followed by F 2 breakdown to relatively low fitness (delayed outbreeding depression); these effects are more extreme in the GAM than in the IFM.

It should also be noted that F 2 breakdown among long-selfed lineages at a purged equilibrium makes them reproductively isolated from each other as they are from their outcrossed ancestors or offspring. This would happen if even a single character conforms to the assumptions of the GAM, with the rest more nearly obeying the IFM.

Predominantly selfing species that exist for long may therefore undergo the progression of reproductive isolation classically observed between young species compared to that between old species, beginning with F 1 heterosis with F 2 breakdown and culminating in F 1 and F 2 sterility and inviability [START_REF] Dobzhansky | Genetics of the Evolutionary Process[END_REF][START_REF] Coyne | Speciation[END_REF].

Evolution of selfing rate near outcrossed and purged equilibria

At selfing rates below the purging threshold, r < r, inbreeding depression is an important genetic constraint which, in combination with ecological constraints, controls evolution of the selfing rate by small genetic steps [START_REF] Lande | The evolution of self-fertilization and inbreeding depression in plants. I. Genetic models[END_REF][START_REF] Johnston | Correlations among fertility components can maintain mixed mating in plants[END_REF][START_REF] Porcher | Evaluating a simple approximation to modeling the joint evolution of self-fertilization and inbreeding depression[END_REF]Devaux et al. 2013). Both the GAM and IFM conform then to the rule that in the absence of ecological constraints selection favors small genetic changes that decrease the selfing rate when the total inbreeding depression exceeds 0.5.

In contrast, in both models inbreeding depression is misleading for understanding evolution of selfing rates above the purging threshold, r > r (or at any purged equilibrium in the GAM). At selfing rates just above the purging threshold the inbreeding depression is reduced or even negative (as in the GAM in Fig 5B ) because most outcrossed survivors are F 1 hybrids between long-selfed lineages (with about half the total genetic variance of their long-selfed parents since f 0 = 0). The progeny of these F 1 outcrosses suffer F 2 breakdown (or delayed outbreeding depression) to low fitness; this strongly selects against outcrossing and favors increased selfing. The models therefore indicate that in a large population the evolution of predominant selfing, with purged quantitative genetic variance, is an irreversible evolutionarily absorbing state [START_REF] Bull | On irreversible evolution[END_REF][START_REF] Charlesworth | Some evolutionary consequences of deleterious mutations[END_REF][START_REF] Galis | Dollo's law and the irreversibility of digit loss in Bachia[END_REF], supporting the view of [START_REF] Stebbins | Self fertilization and population variability in higher plants[END_REF][START_REF] Stebbins | Flowering Plants. Evolution Above the Species Level[END_REF] that predominantly selfing plant species generally occupy terminal branches in plant phylogenies.

Stebbins also inferred that highly selfing species have an increased extinction rate and do not persist long in phylogenetic time. The present theory shows that at a purged equilibrium the long-selfed core of the population has less genetic variance than for a population at the outcrossed equilibrium. This implies that in a constant environment highly selfing populations have a higher mean fitness at equilibrium than under random mating. Extreme environments, such as on the edge of a species range where new selfing species often arise from outcrossers [START_REF] Wright | Evolutionary consequences of selffertilization[END_REF], may exert stronger stabilizing and/or directional selection. Under this scenario, or in a changing environment with high temporal variance and predictability of the optimum phenotype, highly selfing (purged) populations with less genetic variance than at the outcrossing equilibrium will have lower mean fitness through time and thus a higher extinction rate [START_REF] Lande | The role of genetic variability in adaptation and population persistence in a changing environment[END_REF][START_REF] Lande | Stochastic Population Dynamics in Ecology and Conservation[END_REF]. Our model thus provides a theoretical foundation for the classical view of [START_REF] Stebbins | Self fertilization and population variability in higher plants[END_REF][START_REF] Stebbins | Flowering Plants. Evolution Above the Species Level[END_REF] confirmed by recent empirical findings [START_REF] Takebayashi | Is self-fertilization an evolutionary dead end? Revisiting an old hypothesis with genetic theories and a macroevolutionary approach[END_REF][START_REF] Goldberg | Species selection maintains self-incompatibility[END_REF][START_REF] Goldberg | Tempo and mode in plant breeding system evolution[END_REF][START_REF] Wright | Evolutionary consequences of selffertilization[END_REF][START_REF] Igic | Is self-fertilization an evolutionary dead end?[END_REF]) that predominant selfing constitutes an "evolutionary dead end".

This formula can be evaluated from the model only for partial selfing, 0 < r < 1. For random mating and complete selfing, r = 0 and r = 1, the equilibrium inbreeding depression can be obtained by calculating the phenotypic distributions of offspring that would be produced by experimental selfing and outcrossing. The F 1 outcross between two long-selfed parents has a total genetic variance about half that of its parents (because f 0 = 0), hence a higher mean fitness than any selfed age class or any other outcrossed subclass. However, any of their (outcross or selfed) F 2 progeny will have increased genetic variance (due to reduction of negative linkage disequilibrium by recombination), substantially diminishing their mean fitness.

G * = 1 w ∞ τ =1 p τ wτ G τ -γ τ G 2 τ + p 0 ∞ i=0 ∞ j=0 p 0ij w0ij G 0ij -γ 0ij G 2 0ij
where the mean fitness in the population is the same as in eq. (3b). The purged equilibrium (dashed lines) has lower genetic variance, independent of the number of characters, and negative inbreeding depression for r just above the purging threshold. Other parameters as in Fig. 1.

Literature Cited

  rather different genetic variances, creating substantial kurtosis within the outcrossed class, particularly when large negative linkage disequilibrium builds up by selection of different compensatory mutations in long-selfed lineages. The outcrossed class then combines (1) a subclass produced by outcrossing between long-selfed individuals, with about half the total genetic variance of their parents (since f 0 = 0) and (2) subclasses produced by matings with at least one outcrossed parent in which recombination halfs the negative linkage disequilibrium inherited from their parents, increasing the total genetic variance. We found that pooling these subclasses of outcrossed individuals into a single class, ignoring kurtosis, can create artifactual oscillations of genetic variance in the first two age classes. To eliminate this artifact, we analyze selection acting separately on all possible types of outcrosses according to the selfing age of the parents.

  infinitesimal model to encompass mixed mating. Our infinitesimal model involves an infinite number of loci, with no assumption on the distribution of allelic effects within loci. It does, however, assume a Gaussian distribution of breeding values within each cohort of a given selfing age (justified by the Central Limit Theorem as for the classical infinitesimal model under random mating). The accuracy of this assumption is monitored numerically using the kurtosis of breeding values in the population.[START_REF] Fisher | The correlation between relatives on the supposition of Mendelian inheritance[END_REF] infinitesimal model concerns an infinite population with an infinite number of loci each having an infinitesimal effect on a quantitative character. Selection then causes no change in allele frequencies at any locus, although it can change the mean phenotype, and the linkage disequilibrium among loci[START_REF] Bulmer | The effect of selection on genetic variability[END_REF]. The total genic variance, V , in the population thus remains constant, but the total genetic variance, G evolves because selection and recombination change the total linkage disequilibrium variance among loci, C. Fisher's infinitesimal model for an infinite population thus does not require mutation to maintain genetic variance.

  neutral model of partial selfing in an infinite population with no selection or mutation, assuming normality of breeding values within selfing age cohorts, κ = 3( 1-r 2 /4 1-r/4 ) (Appendix 2). Under random mating or complete selfing the breeding value in the population is normal. The maximum kurtosis of 3[1 + (2 -√ 3) 2 ] ≈ 3.215 occurs at selfing rate r = 2(2 -√ 3) ≈ 0.536. Thus in Wright's neutral model of partial selfing in an infinite population the deviation from normality of breeding values at equilibrium must be rather small.

  Genetic variance after selection. With multiple loci and a high selfing rate, the average genetic variance in the population G blows up due to a bimodal distribution of selfing ages, with long-selfed lineages dominating. This blow up is caused by accumulation of linkage disequilibrium in the long-selfed lineages. The large genetic variance expressed by recombination and the decay of linkage disequilibrium in the early selfing age classes greatly reduces their mean fitness. With sufficiently high selfing rates the first few selfing age classes are rare, because of the low outcrossing rate and their greatly reduced fitness. The early selfing ages are largely eliminated by selection, and the genetic variance in the population as a whole after selection is

Figure 1 :

 1 Figure 1: Purging of genetic variance under continued selfing in the GAM for each of 25 identical independent characters under stabilizing selection: components of genetic variance (top) and mean fitness (middle) as functions of selfing age, and distribution of selfing ages (bottom) for population selfing rates below (r = 0.78, left) or above (r = 0.8, right) the purging threshold. Note the different scales in panels A and B. In panel B the genetic variance G at young selfing ages, and genic variance V and covariance C, are not at equilibrium; their values depend on the number of generations simulated (here 500,000). In panels E-F dotted lines represent the distribution of selfing age classes in Wright's neutral model. Other parameters: E = 1, σ 2 m = 0.001, n = 10 and ω 2 = 20.

Figure 2 :

 2 Figure 2: Purging of genetic variance under continued selfing in the IFM for each of 25 identical independent characters under stabilizing selection. A. Total genetic variance before and after selection. Panels and parameters as in Fig. 1 but for the IFM with V (0) = 1. At selfing rates above the purging threshold, r > r, the IFM shows smaller changes in genetic variance and mean fitness as a function of selfing age than in the GAM (Fig. 2. B,D compared to Fig. 1B,D) but still displays similar patterns of selfing age distribution (Fig. 2E,F compared to Fig. 1E,F).

Figure 3 :

 3 Figure 3: Equilibrium genetic variance as a function of selfing rate for each of 25 identical uncorrelated characters under stabilizing selection in the GAM. A. Total genetic variance before and after selection, G and G * . B. Kurtosis of breeding values, κ, before and after selection at the outcrossed equilibrium for r below the purging threshold, at the purged equilibrium for r above the purging threshold, and in Wright's neutral model. Other parameters as in Fig. 1.

Figure 4 :

 4 Figure 4: Equilibrium genetic variance as a function of selfing rate for each of 25 identical uncorrelated characters under stabilizing selection in the IFM. A. Total genetic variance before and after selection. B. Kurtosis in breeding value, κ, before and after selection, and in Wright's neutral model. Other parameters as in Fig. 2.

Figure 5 :

 5 Figure 5: Equilibrium genetic variance after selection, inbreeding depression and population mean fitness, as a function of population selfing rate, for different numbers of characters in the GAM. When the selfing rate is below the purging threshold, two stable equilibria exist. The outcrossed equilibrium has relatively large genetic variance and inbreeding depression nearly independent of selfing rate (solid lines); it is reached when the population has initially low genetic variance and low linkage equilibrium.

Figure 6 :

 6 Figure 6: Equilibrium genetic variance before selection in the IFM, along with inbreeding depression and population mean fitness, as functions of population selfing rate, for different numbers of characters m. Other parameters as in Fig. 2.

Figure 7 :

 7 Figure 7: Equilibrium genetic variance before selection in the IFM, inbreeding depression and population mean fitness, as a function of population selfing rate, for different numbers of characters m with no pleiotropy subject to correlational selection between m/2 independent pairs of characters. In the selection matrix (eq. 2) off-diagonal elements for pairs of characters under correlational selection are 0.5 times the diagonal elements (ω 2 ). The equilibrium genetic covariance between pairs of characters under correlational selection is caused by linkage disequilibrium, and B/G represents their genetic correlation. Other parameters as in Fig. 2.

  Figure1

Acknowledgments

Support was provided by a Royal Society Research Professorship and a grant from the Balzan foundation to R. L., the French CNRS program PICS grant #5273 to E. P., and time on the computing cluster (UMS 2700 OMSI) at MNHN. We thank B. and D.

Charlesworth, W. G. Hill, and J. Pannell for discussions.

Discussion

 [START_REF] Wright | Systems of mating. II. The effects of inbreeding on the genetic composition of a population[END_REF][START_REF] Wright | Evolution and the Genetics of Populations[END_REF]showed that for a neutral model of additive genetic variance, inbreeding increases the equilibrium genetic variance in proportion to the average inbreeding coefficient in the population, so that complete selfing doubles the equilibrium genetic variance in comparison to random mating. [START_REF] Lande | The influence of the mating system on the maintenance of genetic variability in polygenic characters[END_REF] found for regular systems of mating, in which every individual performs the same type of non-random mating, that the equilibrium genetic variance maintained by mutation and stabilizing selection is independent of the mating system. In contrast, [START_REF] Charlesworth | Quantitative genetics in plants: the effect of the breeding system on genetic variability[END_REF] compared maintenance of quantitative genetic variance by mutations under stabilizing or purifying selection, concluding that selfing substantially reduces the equilibrium genetic variance compared to random mating. These disparate results are reconciled in our models of mixed mating, showing that in both the GAM and IFM the equilibrium genetic variance remains nearly the same as under random mating for selfing rates up to the purging threshold, above which a qualitative change occurs. Our asymptotic analysis of the total genetic variance under complete selfing agrees qualitatively with [START_REF] Charlesworth | Quantitative genetics in plants: the effect of the breeding system on genetic variability[END_REF], although their models neglected linkage disequilibrium which plays an important role in purging quantitative variance in both the GAM and IFM.

In the GAM two possible stable equilibria exist for the genetic variance as a function of the population selfing rate, r, with regions of sharply distinct dynamics separated by a purging threshold, r. The outcrossed equilibrium has genetic variance nearly equal to that under random mating, and exists only for selfing rates below the purging threshold. The purged equilibrium has lower genetic variance in the highly selfed stable core of the population, and exists for all selfing rates; in the highly selfed core, the total genetic variance remains constant but the total genetic variance in the population blows up because the genic variance and (negative) linkage disequilibrium in the core increase indefinitely. Initial outcrossing between long-selfed lineages decreases the genetic variance by half (since f 0 = 0), resulting in outcrossed F 1 hybrid vigor ("heterosis"); subsequent selfing or outcrossing of these F 1 allows recombination to express the built-up genic variance previously hidden by negative linkage disequilibrium, producing F 2 breakdown in fitness. Thus, the highly selfed core of a population at a Appendix 1

Truncation of the selfing age distribution

For selfing ages larger than some value L, the genetic parameters within age classes should be nearly equivalent, so that all selfing age classes can be lumped into a single terminal class of age L or larger, which then becomes the upper limit instead of ∞ in numerical summations. Additional recursion formulas are then required for the frequency and genetic composition in the terminal age class.

Recursions for the genic variance and covariance (due to linkage disequilibrium), and the inbreeding coefficient, in the terminal selfing age class are

The recursion for the terminal age class frequency is

where p * τ = wτ p τ / w and w = L τ =0 p τ wτ .

Appendix 2

Population statistics

Kurtosis. By assumption each selfing age cohort, and the population as a whole, has its mean breeding value at the optimum, so there is no skew (asymmetry) in the population. However, because the genetic variance within selfing age cohorts varies among selfing ages, τ , the breeding value in the population will be leptokurtic. The standardized kurtosis of breeding value in the population is the weighted average fourth central moment within cohorts, divided by the square of the population variance in breeding value. For a normal distribution with variance σ 2 the fourth central moment is 3σ 4 , with a standardized kurtosis of κ = 3. Assuming that the distribution of breeding value within each selfing age (and within subclasses of the outcrossed class) is normal, the standardized kurtosis of breeding value in the population is

For comparison, the kurtosis of breeding value in the population can be derived for [START_REF] Wright | Systems of mating. II. The effects of inbreeding on the genetic composition of a population[END_REF][START_REF] Wright | Evolution and the Genetics of Populations[END_REF] classical model of partial selfing with no selection and no mutation. The recursion for the inbreeding coefficient (eq. 13a) with f 0 = 0 has the solution f τ = 1 -(1/2) τ . In the absence of selection the equilibrium selfing age distribution is p τ = (1 -r)r τ . The average inbreeding coefficient in the population is then f = ∞ 0 p τ f τ = r 2-r in agreement with [START_REF] Wright | Systems of mating. II. The effects of inbreeding on the genetic composition of a population[END_REF][START_REF] Wright | Evolution and the Genetics of Populations[END_REF]. With no selection under any amount of outcrossing (r < 1), the population eventually approaches linkage equilibrium, C τ = 0. The equilibrium kurtosis of breeding value in the population is

1-r/4 ). In Wright's neutral model, as in the GAM, at r = 0 or 1 the breeding value in the population is normal (κ = 3). But in the IFM at r = 1 G = 0 (eq. 15) so κ is not defined.

Inbreeding depression. The total inbreeding depression in the population caused by selfing, δ, is one minus the ratio of mean fitness of selfed individuals divided by the mean fitness of outcrossed individuals,


Johanne Bakalara
email: johanne.bakalara@irisa.fr

Thomas Guyet

Olivier Dameron

Emmanuel Oger

André Happe

Temporal models of care sequences for the exploration of medico-administrative data

Keywords: Temporal logics, Description logic, Ontologies, Pharmacoepidemiology, SNDS

Introduction

Pharmacoepidemiology studies uses and effects of drugs on population in real life including benefits and risks. As such, pharmacoepidemiology deals with positive impact (i.e. prevention of disease) as well as safety concerns. Patients with health events are identified, according to their individual characteristics, their treatments and their concomitant treatments.

Collecting information to answer one epidemiological question requires a lot of time and is very expensive.

Medico-administrative databases are a potential alternative which are attractive because of their large population coverage and their availability. For instance, the SNDS1 [START_REF] Tuppin | Value of a national administrative database to guide public decisions: From the système national d'information interrégimes de l'assurance maladie (sniiram) to the système national des données de santé (snds) in france[END_REF] database contains individual information of French patients: age, sex, location; and health reimbursement information: drug deliveries, medical acts or medical visits and hospitalisations (date of arrival, leaving date, diagnosis code) but it does not contain medical reports.

The interest of using this medico-administrative database has been demonstrated by the suspension of benfluorex [START_REF] Weill | Benfluorex and valvular heart disease: a cohort study of a million people with diabetes mellitus[END_REF]. It was the first large-scale pharmacoepidemiology study in France that was possible thanks to information contained in the SNDS.

The study shows that diabetic patients exposed to Benfluorex have an higher risk of hospitalisation for heart disease than the unexposed diabetics in the following years. Epidemiologists use the SNDS database to find patients that experienced some medical events of interest. It is worth noticing that information has been collected for administrative purposes (care reimbursements), but not for medical ones. As a consequence, the medical content associated with medical events is relatively poor. Let us consider a patient having a Venous Thromboembolism (VTE). The database does not contain the information that the patient has the disease VTE. This information has to be deduced from the information about its care deliveries. For instance, a consultation with a vascular specialist and deliveries of anticoagulant drugs. These pieces of information are available in the data in a structured format. Compared to medical record in hospital, the SNDS database has less medical information, but it does not required sensitive text analysis.

All medical events that are stored in the database compose the care trajectory of a patient. Figure 1 illustrates a care trajectory. The challenge of the epidemiologist is to define selection criteria that would reconcile those actual patient information with the medical semantic. The definition of these criteria composes a health pattern called care pathway illustrated Figure 2. A possible care pathway of patients with VTE2 is a consultation with a vascular specialist or at the hospital with an Doppler imagery act (CCAM3 EDQM001 code) followed by at least 3 deliveries of anticoagulant drugs within 4 months (with the first delivery of anticoagulant in the week after consultation/Doppler). The care trajectory denotes the actual medical events stored in the database while a care pathway is a pattern of medical events that describe the pathology profile as illustrated on Figures respectively 1 and 2.

We aim at enabling epidemiologists to query the database of care trajectories with such kind of complex descriptions of care pathway. The complexity is twofold:

• use of ontological concepts (Doppler imagery act/anticoagulant/vascular specialist):

the code of the medical act is given, but the code for anticoagulant drugs is not precisely given. Anticoagulants refer to a class of drugs that is described in the ATC taxonomy.4

• use of temporal constraints (in the week after/within 4 months): the temporal order of cares, numerical duration/delays specifies the temporal organisation of the events.

Our work focuses on having expressive temporal constraints in the description of care pathways. Temporal constraints enable to discriminate care trajectories of interest from care trajectories with same events but presenting different delays (e.g. an Doppler imagery act followed by an anticoagulant delivered several weeks later does not witness a VTE). This may help to specify care pathways that match only the desired care trajectories. It is worth noting that drugs deliveries, medical acts and hospital stays are available in the database with timestamps.

The challenge we face is summed up in Figure 3. Care pathways are on the left side and care trajectories are on the right side. The overall objective is to bridge the two sides. On the top in red, epidemiologists formulates a medical hypothesis. In blue, at the second level, epidemiologist describes care pathways with the available health information and a database stores health information. In green, at the third level, data and queries are specified in a formalism. This latter has to express a maximum of information to represent the complexity of the data (care trajectories) and of the care pathway. Green and blue levels are intertwined. Today, the tools for querying SNDS (third level) lack of expressiveness and constraint epidemiologists for their specification (SQL, SAS, etc...). In this article, we show that the problem of querying patient care trajectories with temporal medical patterns (the care pathways) may be addressed, sometimes only partially, by tools or formalisms coming from different computer science fields: complex event processing, formal reasoning, knowledge representation and database. In this preliminary work, our objective is to identify strengths and weaknesses of these

In the next sections, we states formally of problem of querying a set of care pathways, then Section 3 reviews related work. The last section illustrates some of approaches of the related work on a case study.

Problem statement

The objective of this work is to propose a formal framework that would design a well-founded and efficient tool for querying care trajectories in the context of pharmacoepidemiology.

Generally speaking, let T = (T i ) i∈[n] be a set of n care trajectories and ϕ a care pathway abstract description. ϕ holds in a care trajectory T ∈ T , denoted T |= ϕ, if and only if the care trajectory contains the care pathway. The formalisation problem is threefold:

• define a formalism to model care trajectories, T , which represent the SDNS data

• define a formalism to model care pathways, ϕ, which specifies an abstract care pathways

• define a computational model that can evaluate whether T entails ϕ: T |= ϕ.

As we noticed in the introduction, specifying care pathways requires to manipulate: temporal concepts (time constraints and time window), medical concepts and knowledge (ontologies). The ideal formal framework should capture these dimensions, enable intuitive queries to be expressed for a wide range of pharmacoepidemiological studies and be computationally efficient.

It is of paramount importance to base choices on solid theoretical foundations. Expressiveness and efficiency are known to be antagonist objectives [START_REF] Levesque | Knowledge representation and reasoning[END_REF]. A well-founded approach would be the basis for proposing long-term solutions, make possible future improvements and facilitate its application to a broad range of contexts (i.e., various databases, queries). This section presents four families of formalisms to answer the problem: model checking; Complex Event Processing (CEP); temporal databases; and Knowledge Representation and Reasoning (KR). The work that we mentioned addresses the three problems mentioned in Section 2 at the same time. The formalism should represent data (care trajectories), query (care pathway) and to compute the answers of the queries on data. The last two families have been more explored in medical context [START_REF] Combi | Temporal information systems in medicine[END_REF] than the two others.

Model checking [START_REF] Clarke | Model checking[END_REF] verifies if a model satisfies a property or a formula. This research line is interested in representing dynamic systems with formal temporal formalisms (discrete event models, M, describing how the system evolves) to prove some properties specified by formula ϕ. A formula ϕ is true if and only if ϕ is true for any traces5 that can be generated from the model M. The most common formalisms for formula in Model Checking are LTL (Linear Temporal Logic), CTL (computation tree logic) or MTL (metric temporal lgic) which is the temporal extension of LTL.

To apply such methods in our context, the events of care trajectory are represented by one finite trace of the system (and there is no system model in our case) and the care pathways is represented by a temporal logic formula. The care trajectory is selected if and only if the trace satisfies the formula. These methods are interesting because they provide formal results (expressiveness, completeness, equivalences) on the representation of timed systems, but they don't manage neither reasoning nor ontological representation. In related medical domain, model checking has been used to study the compliance of care pathways [START_REF] Bottrighi | Adopting model checking techniques for clinical guidelines verification[END_REF].

Complex Event Processing (CEP) [START_REF] Giatrakos | Complex event recognition in the big data era[END_REF]) is a research line that aims at processing log-streams with patterns. Log-streams are streams or sequences of timed events. The CEP processes these logs to detect or to locate complex events (or patterns) defined by the user. This domain defines formalisms that aim at being very efficient to process streams and expressive to specify patterns. Temporal constraint networks [START_REF] Cabalar | Temporal constraint networks in action[END_REF] or Chronicles [START_REF] Dousson | Chronicle recognition improvement using temporal focusing and hierarchization[END_REF] are simple temporal models that are interesting for their graphical representation, but are limited to simple relational events. Some more complex formalisms, e.g. ETALIS [START_REF] Anicic | ETALIS: Rule-based reasoning in event processing[END_REF] or logic-based event recognition [START_REF] Giatrakos | Complex event recognition in the big data era[END_REF], propose very expressive representations of complex events, including reasoning techniques (including ontologies) which enrich the capabilities of CEP.

In our context, care trajectories are logs, and care pathways are the complex events. We are not interested in the stream dimension of these approaches, but their formalisms to represent complex events may be adapted in the context of static logs.

Temporal databases (Snodgrass, 1992) extend the notion of database to timestamped data. Databases issued data representation problems but also specific querying language problems. We gather in this family the temporal extension of relational databases (e.g. TSQL) but also web semantic approaches which combine query languages (e.g. SPARQL) and expressive description languages. Care trajectories are facts in the temporal database and the querying of care pathways becomes a problem of specifying care pathways in the query language. [START_REF] Rivault | queryMed: Semantic web functions for linking pharmacological and medical knowledge to data[END_REF] shown that semantic web is an interesting approach for our problem, but does not explicitly address the problem of timed queries.

Finally, Knowledge Representation and Reasoning (KR) [START_REF] Levesque | Knowledge representation and reasoning[END_REF] is "the study of how what we know can at the same time be represented as comprehensibly as possible and reasoned with as effectively as possibly". In this research domain, temporal KR is focused on representing and reasoning about time. It gives rise to several logics [START_REF] Long | A review of temporal logics[END_REF], for instance: Allen's logic, McDermott's logic, Event Calculus or Halpern & Shoham's logic.

KR is a general framework to study how to represent care trajectories and how to model reasoning-based queries on care pathways. Approaches from the other families may be represented with appropriate logics. Studying KR formalisms seems of paramount importance as it provides common foundations to compare various approaches. Description Logic (DL) [START_REF] Baader | The description logic handbook: Theory, implementation and applications[END_REF] is a KR formalism allowing ontology-mediated query answering (OMQA) [START_REF] Bienvenu | Ontology-mediated query answering: harnessing knowledge to get more from data[END_REF]. Artale et al. [START_REF] Artale | Ontology-mediated query answering over temporal data: A survey[END_REF] present a temporal extension of DL that may be suitable for our problem. For instance, [START_REF] Shankar R | Knowledge-data integration for temporal reasoning in a clinical trial system[END_REF] developed a tool based on OWL for research data management with a temporal reasoning in a clinical trial system.

4 Comparison of approach on a use case

Rational

This section introduces a real use case. In this example, pharmacoepidemiologists want to select patients with Venous Thromboembolism (VTE) from the data contained in the SNDS. Venous thromboembolism, i.e. deep vein thrombosis (DVT) or pulmonary embolism (PE), is a frequent and potentially fatal disease [START_REF] Oger | Incidence of venous thromboembolism: a community-based study in western france[END_REF][START_REF] Delluc | Real-world incidence of cancer following a first unprovoked venous thrombosis: Results from the EPIGETBO study[END_REF]. This requires to survey how many people are concerned, if the number of patients increased and if a specific drug has an impact. The difficulty for epidemiologists lies in the description of the care pathways that will accurately identify VTE from the SNDS data. The description below describes two care pathways that physicians proposed to identify VTE (referring to Figure 3, this description is between the red and blue levels).

In clinical practice facing a clinical suspicion of VTE, physicians first prescribe anticoagulant and then confirm or not the diagnosis through specific medical acts: for instance Doppler ultrasonography or CT scan. Patients with suspected PE are often hospitalized whereas patients with suspected DVT are managed on a ambulatory basis. If the suspicion is confirmed, anticoagulant deliveries continues for 3 to 12 months or sometimes longer duration. Hence, diagnosis (through medical act) is preceded or followed by anticoagulant initiation within a time window of at most 0 to 2 days, keeping in mind that PE suspicion leads to hospitalisation during which medical act to confirm the diagnosis are performed and then anticoagulant is observed only after the patient comes back home.

Through these observations, pharmacoepidemiologists identified the following two care pathways to detect patients with VTE form SNDS data (referring to Figure 3 this description is in the blue level):

1. A diagnosis (DVT or PE) or a medical act (Doppler or CT scan) during or prior to anticoagulant(AC) deliveries for 1 to 2 days and delivery lasts a minimum of 3 months and a maximum of 12 months (sometimes longer). Each delivery is separated by 0 to 2 months.

2. A diagnosis PE during an hospitalisation followed by AC delivery.

These care pathways contain sequential order but also time constraints between events (for instance number of days) or duration of events (time window). Searching for such patterns requires high expressiveness that make databases query languages (SQL, TSQL) practically difficult to use. Next section illustrates alternative solutions to represent the green level (see Figure 3).

This use case illustrated the problem of formalizing care pathways of patients suffering from VTE. But, for sake of generality, our formalism has to specify the care pathway patterns of as many case studies as possible.

IA & Santé 2019

Comparison of different models

This section aims at illustrating the comparison of four formalisms, Description Logic (DL), Chronicles, LTL and MTL to discuss about their power of expressivity in our case of study.

DL supplies ontologies query answering technologies for Ontology-Based Data Access (OBDA) and is well known to represent medical data. Indeed, Description Logic is used to describe and reason about concepts on data. Reasoning with Description Logic is performed in three steps. The first one consists in defining the data and the data form (called ABox) which contains knowledge at the instance level: a set of assertion defining concepts, roles and a countably infinite set of individuals names. Concepts with individuals names and roles with individuals names are forming atoms.

The second step consists in define a base of knowledge (TBox) which is a set of concepts inclusions. Concept inclusions represent a hierarchy of concepts. For instance, the concept B01AF 01 representing anticoagulant drugs is the leaf concept in the hierarchy of concept modeling the ATC classification:

• B: Blood and blood forming organs

• B01: Antithrombotic agents • B01AF: Direct Xa inhibitor • B01AF01: Rivaroxaban And considering the CCAM (medical acts) code for the Doppler: EDQM001 (iliac and lower limb arteries) we could construct the following ABox, where Pierre and Paul are patients and n i are dates of medical events:

B01(Patient1, t 1 ) B01(Patient1, t 2 ) P E(Patient1, t 3 ) B01AF 01(Patient2, t 1 ) B01AF 01(Patient2, t 3 ) EDQM 001(Patient2, t 4 ) t 1 < t 2 < t 3 < t 4
And TBox issued from the previous ATC classification:

B01AF 01 B01AF B01AF B01 B01 B
B01AF 01 is the subclass of B01AF which is the subclass of B01 which is the subclass of B. Operators linking concepts are defined depending on the class of DL chosen. For example, the ALC family offers the concept constructors: negation, conjunction, disjunction, existential restriction ∃ and universal restriction ∀.

The third step is to define a query to extract information from knowledge contained in ABox and TBox. Usually, queries are expressed with a first order logic. [START_REF] Artale | Ontology-mediated query answering over temporal data: A survey[END_REF] designed a temporal DL: TQL that extends the standard ontology language: OWL 2QL. It offers the capability of having time as individuals names and to compare whether an atom occurs before another. In our case study, we propose a query to find patients designed by the use case defined in Section 4.1.

∃patientID, t 1 (t ref < t 1 ) ∧ EDQM 001(patientID, t ref ) ∧ B01(patientID, t 1 ) ∧t 2 (t 2 < t 1 ) ∧ B01(patientID, t 2 ) ∧t 3 (t 3 < t 1 ) ∧ B01(patientID, t 3 )
We find patients having three deliveries of B01 (and subclasses of B01) in a row and the first one is delivered after an Doppler. The temporal information that can be represented here is limited to the order of the events. Delay or duration can not be specified.

To express queries with temporal delays, the formalism of chronicle represents a care pathway as a temporal constraint graph. Chronicles allow the expression of sequential order of events with temporal constraints such as interval of time. Furthermore, negative time in the interval expresses that an event may occur before or after another one. Figure 4 specifies patients having at least three anticoagulant deliveries separated by 0 to 60 days, and a diagnosis DVT before, after or during deliveries. DVT occurs 92 days earlier or one day after the third delivery. However, we can not explicitly restrict the number of deliveries to 12 months as defined in the use case. We also cannot use the notion of no event (event does not occur). Model checking offer the possibility to express no event and can be used as queries. Such as an example we propose the following LTL formula as an example applied to our case of study:

(♦D B01 ∧ (♦D B01 ) ∧ (♦D B01 ) ∧ (♦(DV T ∨ P E))
The LTL formula represents a care pathway with at least three deliveries and a diagnosis DVT or PE. We literally read it: in the future (♦), there is the delivery of B01 and (∧) it is followed ( ), in the future, by the delivery of B01 and it is followed, in the future, by a delivery of B01 and, in the future, there are the diagnosis DVT or (∨) PE). LTL only contains order between events and doesn't contains time constraints. It is quite limited for our problem, so we refer to its temporal extension MTL. The MTL formula adds the capability to express quantitative temporal constraints. We propose the following MTL formula as an example applied to our case of study:

♦(DV T ∨ P E) ⇒ ((♦ [0 2] D BO1 ) ∧ (♦ [0 60] D BO1 ) ∧ (♦ [0 60] D BO1 ) ∧ (♦ ≥365 ( (¬D BO1 )))
It represents a care pathway with a DVT or PE followed between 0 to 2 days after by three AC deliveries separated between 0 to 60 days, and no deliveries occur after 365 days. MTL can explicitly restrict the number of deliveries and temporal constraints but the notion of sequences is manually found by the multiple use of ♦. For deep understanding of notations, we refer the reader to [START_REF] Bouyer | On the expressiveness of tptl and mtl[END_REF].

From computational point of view, Chronicles may be very space/time-efficiently to be recognizes in care trajectories. Simple LTL formula would also be space/time-efficient to check but it is expressively poor. In contrast, MTL is known to be undecidable. It is a theoretical limitation but, not necessary a practical constraint [START_REF] Ouaknine | On the decidability of metric temporal logic[END_REF].

Conclusion

In this article, we introduced the context of pharmaco-epidemiological studies with medicoadministrative databases and the challenge to query such databases with medical questions. It consists in finding patients satisfying a medical pattern that we want to be expressive.

To express medical patterns, we use the formalism of Description Logic to describe data and include medical knowledge issued from the available classification (e.g. ATC, CCAM). To query these data, we compared solutions issued from first order logic, Chronicles and MTL. It shows that none of them is enough expressive to correctly translate all the desirable constraints of the care pathway. The future work consists of extending these formalisms and to study efficiency issues by testing them through many cases of study on the SNDS.

Figure 2 :

 2 Figure 1: Example of a care trajectory.

Figure 3 :

 3 Figure 3: From medical study issue to the care trajectory querying

  Figure 4: Chronicles

SNDS: French National System for Health Data (previously SNIIRAM). The SNDS is the world largest medico-administrative database with a population coverage close to 99%.

This description is an illustration for pedagogical purposes, we refer to the use case for a medical description of VTE.

CCAM: Classification Commune des Actes Médicaux

ATC: Anatomical Therapeutic Chemical Classification System.

Sequences taht register the set of atomic propositions that are valid along the execution, cf. part 3.2.2 in[START_REF] Baier | Principles of model checking[END_REF]