T Preece

Y Mao

The evolution of gynodioecy on a lattice

Keywords: sexual reproduction, lattice, Gynodioecy

Gynodioecy is a breeding system in plants where populations consist of hermaphrodites and females. The females result from a genetic mutation which impairs pollen production in hermaphrodite plants. Most previous models for the evolution of gynodioecy do not take into account any spatial detail, which might be expected to play an important role in populations with short range interactions caused by poor or no locomotion.

In this article we present a generalised mean-field analysis (which ignores any spatial detail), together with stochastic spatial simulations, to investigate the spatial effect on the evolution of gynodioecy. We show that, in a population of hermaphrodites where male sterility is caused by a dominant allele in a nuclear gene, mean-field calculations greatly underestimate the reproductive advantage females require to become viable under spatial constraints. This suggests that gynodioecy is less likely to evolve in plants with more localised pollination and seed setting. This may have implications for the evolution of dioecy, a breeding system in plants where the population consists of males

Introduction

Gynodioecy is a breeding system in plants where populations consist of hermaphrodite and female individuals. It is both a common and widespread polymorphism describing approximately 7% of all flowering plants [START_REF] Richards | Plant breeding systems[END_REF]. Some examples of species that exhibit gynodioecy are the wild strawberry, Fragaria virginiana, [START_REF] Ashman | Determinants of sex allocation in a gynodioecious wild strawberry: implications for the evolution of dioecy and sexual dimorphism[END_REF] and Cucurbita foetidissima (Cucurbitaceae) [START_REF] Kohn | Sex-ratio, seed production, biomass allocation, and the cost of male function in cucurbita-foetidissima (cucurbitaceae)[END_REF]. The evolution of gynodioecy also has implications for the evolution of dioecy, as gynodioecy is thought to be a route to dioecy.

A well-known dioecious species is American holly, Ilex aquifolium (Aquifoliaceae) [START_REF] Obeso | Sex ratios, size distributions, and sexual dimorphism in the dioecious tree ilex aquifolium (aquifoliaceae)[END_REF].

Gynodioecy occurs as a result of a genetic mutation which impairs pollen production in hermaphrodite plants. This makes gynodioecy an excellent system for studying the interplay of genetic architecture and ecology in evolution [START_REF] Bailey | A field guide to models of sex-ratio evolution in gynodioecious species[END_REF]. Genetic mutations can impair pollen production in several ways but still allow normal female reproduction [START_REF] Chaudhury | Nuclear genes-controlling male-fertility[END_REF].

The genetic basis for male sterility can significantly affect whether a mutation for male sterility can become established. Nuclear genes are inherited through both parents whereas cytoplasmic genes are only inherited from the mother. It has been shown that females would need to be at least twice as fertile as hermaphrodites (produce twice as many offspring) for a nuclear male sterility gene to become established [START_REF] Lewis | Male sterility in natural populations of hermaphrodite plants[END_REF]. However, cytoplasmic male sterility can become established when the fertility of the female is only slightly larger than that of the hermaphrodite [START_REF] Lewis | Male sterility in natural populations of hermaphrodite plants[END_REF]. This is because a cytoplasmic gene is not inherited through male function hence has nothing to lose if pollen production is stopped. In some species, male sterility is complicated further by the evolution of restorer genes which suppress the deleterious effects of cytoplasmic male sterility [START_REF] Schnable | The molecular basis of cytoplasmic male sterility and fertility restoration[END_REF].

Most previous models of inherited male sterility assumed that female fitness is dependent on the availability of pollen in the environment but hermaphrodite fitness is not [START_REF] Lewis | Male sterility in natural populations of hermaphrodite plants[END_REF][START_REF] Lloyd | Evolutionarily stable sex-ratios and sex allocations[END_REF], i. e. hermaphrodites are self-compatible. However, incompatibility recognition systems preventing self-fertilization have evolved several times in independent lineages of Angiosperm plants [START_REF] Charlesworth | Plant selfincompatibility systems: a molecular evolutionary perspective[END_REF]) and many plants have been observed as self-incompatible [START_REF] Igic | The distribution of plant mating systems: Study bias against obligately outcrossing species[END_REF]. If the case of a selfincompatible species is considered, cytoplasmic inherited male sterility is only stable if hermaphrodite and female fertility are exactly equal [START_REF] Charnov | The Theory of Sex Allocation[END_REF]. In this unlikely case it would be difficult for the cytoplasmic mutation to invade from small frequencies as it is neutral in terms of selection. For the case of cytoplasmic male sterility in out-crossing hermaphrodites, it was shown that if female fertility is larger than hermaphrodite fertility, females will spread until scarcity of pollen causes the population to become extinct [START_REF] Stewart-Cox | Space mediates co[END_REF].

All of the analyses discussed thus far ignore spatial detail, which becomes important in populations with short range interactions between constituents and limited mixing caused by poor/no locomotion. These conditions would be expected to be particularly prevalent in plant populations. A powerful approach for modelling spatial detail is the use of stochastic spatial simulations, which have been extensively applied in chemical, ecological and sociological systems [START_REF] Durrett | Stochastic spatial models[END_REF][START_REF] Durrett | Stochastic spatial models a users guide to ecological applications[END_REF].

A spatial model of cytoplasmic male sterility in self-incompatible hermaphrodites [START_REF] Stewart-Cox | Space mediates co[END_REF], showed behaviour vastly different from the predictions of non-spatial models. Where non-spatial models predicted that invasion of females would lead to extinction of the population (female fertility > hermaphrodite fertility), the spatial model [START_REF] Stewart-Cox | Space mediates co[END_REF] instead displayed stable nodes, foci, limit cycles or extinction depending on the relative fertility of females to hermaphrodites. This demonstrated that cytoplasmic male sterility can evolve in self-incompatible hermaphrodite populations.

In this article we present a generalised model for gynodioecy that accommodates male sterility conferred by a cytoplasmic gene or a dominant allele in a nuclear gene. This contrasts the aforementioned spatial investigation [START_REF] Stewart-Cox | Space mediates co[END_REF], which solely considered cytoplasmic male sterility. We present a generalised mean-field analysis (which ignores any spatial detail), together with stochastic spatial simulations, to investigate the spatial effect on the evolution of gynodioecy in self-incompatible hermaphrodites.

For established gynodioecious populations we compare the equilibrium frequencies of females in spatial and non-spatial models over a range of parameter values. A knowledge of the equilibrium frequency of females would be of interest when comparing the results of breeding experiments with observations of female frequency in the wild.

In the following section a stochastic spatial simulation is described that is used to simulate the evolution of gynodioecy. A mean-field analysis for the evolution of gynodioecy is presented in section 3 and results from the stochastic spatial simulations are presented in section 4. The results presented for the case of cytoplasmic male sterility show good agreement with previous spatial simulations [START_REF] Stewart-Cox | Space mediates co[END_REF].

Model

The stochastic spatial simulation described here is continuous in time and extends a previous model for the sexual reproductive process [START_REF] Preece | Sustainability of dioecious and hermaphrodite populations on a lattice[END_REF]. The population resides on a two-dimensional square lattice of side length L, with periodic boundary conditions. Each lattice site can be occupied by either a self-incompatible hermaphrodite, a female or be empty.

An occupied site becomes vacant at a rate λ, the death rate. The death rate is considered to be equal for hermaphrodites and females for the purposes of our investigation.

A hermaphrodite produces an offspring at a rate equal to jκ 1 /4, where j is the number of nearest neighbour sites occupied by hermaphrodites and κ 1 is the birth rate. This positive relationship between individual fitness and the density of conspecifics demonstrates the well-known Allee effect [START_REF] Odum | A note on the stable point of populations showing both intraspecific cooperation and disoperation[END_REF]. A nearest neighbour site is then chosen randomly. If the chosen site is empty the offspring is placed on the site. If the site is occupied the offspring is deleted. Similarly, a female produces a hermaphrodite offspring at a rate, jκ 2 /4, and a female offspring at a rate, jκ 3 /4. Thus, a female produces offspring at a total rate, j(κ 2 + κ 3 )/4. Genetics can be explicitly defined in this model for the case of a dominant allele in a nuclear gene that causes male sterility, or a cytoplasmic gene that causes male sterility. Consider a dominant nuclear allele for male sterility, B. Thus, Bb = female, and, bb = hermaphrodite. Since females produce no pollen, no homozygous BB individuals are formed and there are only two genotypes in the population [START_REF] Lewis | Male sterility in natural populations of hermaphrodite plants[END_REF]. κ 2 = κ 3 , represents the case of differential mortality between male and female offspring due to their own sex genotypes [START_REF] Lloyd | Theoretical sex-ratios of dioecious and gynodioecious angiosperms[END_REF].

The case of cytoplasmic male sterility is modelled by setting κ 2 = 0.

Hence, females always produce females and the hermaphrodite population is not increased through mating with females.

Mean-field Analysis

The mean-field kinetics of the hermaphrodite and female populations are described by,

dh dt = h -λ + 3 (hκ 1 + fκ 2 ) (1 -h -f ) 4 (1) df dt = f -λ + 3hκ 3 (1 -h -f ) 4 (2)
where h and f are the densities of hermaphrodites and females respectively.

The derivation of equations ( 1) and ( 2) is explained in appendix A. By scaling time (τ = λt), the system can be simplified by introducing the parameters;

K 1 = 3κ 1 /4λ, K 2 = 3κ 2 /4λ and K 3 = 3κ 3 /4λ.
In the following mean-field analysis, the case of nuclear male sterility (K 1 > 0, K 2 > 0 and K 3 > 0) will be considered first and the case of cytoplasmic male sterility (K 1 > 0, K 2 = 0 and K 3 > 0) will be considered second.

3.1. Nuclear male sterility.

Equilibria for equations ( 1) and ( 2) were found by setting the time derivative equal to zero, then solving the resulting two equations simultaneously.

The nature of these equilibria was characterised by linear stability analysis of the system of equations ( 1) and ( 2). As h ≥ 0, f ≥ 0 and h + f ≤ 1, the region of interest on the (h,f )-plane is the triangle with vertices (h = 0, f = 0), 1) and ( 2)

(h = 1, f = 0) and (h = 0, f = 1). The trivial equilibrium (h = 0, f = 0) is stable provided λ > 0. For (K 1 > 0, K 2 > 0, K 3 > 0), equations (
have four equilibria besides the trivial one. Two equilibria are located on the boundary (f = 0) with,

h 2 = 1 2 1 - K 1 -4 K 1 , (3)
h 3 = 1 2 1 + K 1 -4 K 1 . (4)
The remaining two equilibria are,

h 4 = K 2 K 3 - √ A 2K 3 (K 2 + K 3 -K 1 ) , (5)
f 4 = (K 3 -K 1 ) K 2 K 3 - √ A 2K 2 K 3 (K 2 + K 3 -K 1 ) , (6)
h 5 = K 2 K 3 + √ A 2K 3 (K 2 + K 3 -K 1 ) , (7)
f 5 = (K 3 -K 1 ) K 2 K 3 + √ A 2K 2 K 3 (K 2 + K 3 -K 1 ) , (8)
where

A = K 2 K 3 [4 (K 1 -K 2 -K 3 ) + K 2 K 3 ].
Eigenvalues and eigenvectors from the linear stability analysis of these equilibria can be found in appendix B.

Equilibria (h 2 , 0) and (h 3 , 0) are real and positive only if K 1 > 4. The eigenvalues from the linear stability analysis of these equilibria are only de-

pendent on K 1 and K 3 (see appendix B). If K 1 > K 3 , equilibria (h 2 , 0) and
(h 3 , 0) are a saddle point and attractor respectively. This marks a regime where hermaphroditism is an ESS. Thus, a small number of mutant females cannot invade a hermaphrodite population. Later it will be shown that when

K 1 > K 3
, only equilibria (0, 0), (h 2 , 0) and (h 3 , 0) exist in the region of interest. Typical mean-field dynamics for the regime, where K 1 > 4 and K 1 > K 3 , are shown in figure 1. The unstable manifold of the saddle point in figure 1 lies along the boundary (f = 0). Thus, a population of just hermaphrodites evolves to the trivial equilibrium if the initial density of the population is below the saddle point and evolves to the stable node if the initial density of the population is above the saddle point.

If K 3 > K 1 , equilibria (h 2 , 0) and (h 3 , 0) are an unstable node and saddle point respectively. The unstable manifold of the saddle point at (h 3 , 0) points in to the interior, (h > 0, f > 0). In the absence of females the dynamics of the hermaphrodite population is the same as for the case, K 1 > K 3 . However, a small number of mutant females can now invade the hermaphrodite population. Thus, Hermaphroditism is no longer an ESS.

By inspecting equations ( 5)-( 8) it can be seen that equilibria (h 4 , f 4 ) and (h 5 , f 5 ) are real only if A ≥ 0. Using this condition, it can be shown that equilibria (h 4 , f 4 ) and (h 5 , f 5 ) are positive, real, and distinct if,

K 3 > 4, K 3 > K 1 ,
and

K 2 > S, (9)
where,

S = 4K 3 -4K 1 K 3 -4 . ( 10
)
When they exist in the region of interest, equilibria (h 4 , f 4 ) and (h 5 , f 5 ) are a saddle point and stable node respectively. The condition, in equation ( 9), that K 3 > K 1 , means that the boundary equilibria (h 2 , 0) and (h 3 , 0) are an unstable node and a saddle point respectively. Typical mean-field dynamics for this regime are shown in figure (2). In the absence of females, the dynamics of hermaphrodite population is similar to the case where K 1 > K 3 : a population of just hermaphrodites evolves to the trivial equilibrium if the initial density of hermaphrodites is below the unstable node and evolves to the boundary saddle point if the initial density of hermaphrodites is above the unstable node. However, females can now invade a population of hermaphrodites. The population can then evolve to the stable node (h 5 , f 5 ), which represents a stable coexisting population of hermaphrodites and females. As can be seen in figure(2), mixed populations with low densities evolve to the trivial equilibrium.

The conditions in equations ( 9) also reveal that a coexisting population of females and hermaphrodites can be stable even when a population consisting solely of hermaphrodites is not stable (when K 1 < 4).

Whilst obeying equations (9), reducing K 3 (or increasing K 1 ) brings equilibria (h 4 , f 4 ) and (h 5 , f 5 ) closer to the boundary (f = 0). When K 3 = K 1 , equilibria (h 4 , f 4 ) and (h 5 , f 5 ) collide with the boundary equilibria, (h 2 , 0) and (h 3 , 0) respectively, here a transcritical bifurcation occurs. As discussed above, when K 1 > K 3 , the boundary equilibria (h 2 , 0) and (h 3 , 0) are a saddle point and attractor respectively. Thus hermaphroditism is the ESS.

Decreasing K 2 , whilst obeying the conditions in equation 9, brings equilibria (h 4 , f 4 ) and (h 5 , f 5 ) closer together. At K 2 = S, equilibria (h 4 , f 4 ) and

(h 5 , f 5 ) collide and annihilate. This is a saddle-node bifurcation point. If K 3 > K 1 and K 2 < S; the boundary equilibria, (h 2 , 0) and (h 3 , 0), are a node and saddle point respectively and no equilibrium for coexistence exists.

Under these conditions females could invade the hermaphrodite population.

As the only stable equilibrium is the trivial equilibrium the population will then converge on the trivial equilibrium, resulting in extinction. Typical mean-field dynamics for this regime are shown in figure 3.

All possible evolutionary outcomes for the case K 1 > 0, K 2 > 0 and K 3 > 0 are summarised in table (1). At equilibrium the frequency of females in the population is given by,

f 5 f 5 + h 5 = K 3 -K 1 K 3 + K 2 -K 1 . (11)
Note that equation ( 11) is only valid if the (h 5 , f 5 ) is real and positive, therefore if the conditions in equation 9 are satisfied. Equation 11 differs from the result of [START_REF] Lewis | Male sterility in natural populations of hermaphrodite plants[END_REF] as Lewis's model considers self-compatible hermaphrodites.

K 1 > 4 K 1 < 4 K 1 > K 3 Hermaphroditism Extinction is ESS (only trivial equilibrium) (K 3 > K 1 )

Female invasion Extinction

& (S > K 2 ) leads to extinction (only trivial equilibrium)

(K 3 > K 1 ) & (K 2 > S) Gynodioecy Gynodioecy & (K 3 > 4)
Table 1: Evolutionary outcomes in mean-field analysis, for nuclear male sterility (K 1 > 0, K 2 > 0, K 3 > 0), where K 1 , K 2 and K 3 are the scaled reproductive rates for the following processes: H + H → H, H + F → H, and H + F → F respectively (H represents a hermaphrodite and F a female).

Cytoplasmic male sterility

For the case of cytoplasmic male sterility, the behaviour of equilibria (h 2 , 0) and (h 3 , 0) is the same as in the previous section: for K 1 > K 3 equilibria (h 2 , 0) and (h 3 , 0) are a a saddle point and stable node respectively, for K 3 < K 1 equilibria (h 2 , 0) and (h 3 , 0) are an unstable node and saddle point respectively.

If K 2 = 0 and K 1 > K 3 , hermaphroditism is an ESS and similar dynamics to those shown in figure 1 are observed. For K 2 = 0 and K 3 > K 1 , females can invade a hermaphrodite population but the population then converges on the trivial equilibria and becomes extinct. The dynamics observed for this case are similar to those shown in figure 3.

For the special case K 2 = 0 and K 3 = K 1 , there is a curve of equilibrium points given by,

f = K 1 h -K 1 h 2 -1 K 1 h . ( 12
)
Eigenvalues and eigenvectors from the linear stability analysis of this equilibrium curve are given in appendix C. The first eigenvalue for points on this equilibrium curve is zero with corresponding eigenvector along the curve.

The second eigenvalue is negative for h > 2/K 1 and positive for h < 2/K 1 .

Typical mean-field dynamics for this regime is shown in figure 4. The behaviour observed here for the case K 2 = 0 is in agreement with previous results from similar evolutionary models [START_REF] Stewart-Cox | Space mediates co[END_REF][START_REF] Charnov | The Theory of Sex Allocation[END_REF].

Results of lattice simulation

All simulations were run on lattices of side length L = 100 for 100 generations, where the time in generations is given by τ = λt, and t is in Monte Carlo steps. For hermaphrodite populations (with no females), previous investigations for this lattice model [START_REF] Preece | Sustainability of dioecious and hermaphrodite populations on a lattice[END_REF] have shown that a discontinuous phase transition occurs between an active phase (finite population) and an absorbing phase (extinction) at K 1 ≈ 6. Thus values of K 1 > 6 are required for a sustainable hermaphrodite population.

The simulation results for the case of cytoplasmic male sterility (K 2 = 0)

show good qualitative agreement with previous spatial simulations [START_REF] Stewart-Cox | Space mediates co[END_REF]. Simulations (for K 2 = 0) displayed coexistence of females and hermaphrodites (not shown) for values of K 3 roughly 1.5 times lager than K 1 . Figure 5 shows how behaviour is increasingly cyclic with increased female reproductive advantage (increased K 3 /K 1 ). This contrasts with mean-field kinetics, which predict extinction of the population following the invasion of females. Large vales of K 3 (K 3 > 8K 1 ) result in extinction (not shown)

showing similar behaviour to mean-field predictions.

Sample paths from the lattice simulation for the case of nuclear male sterility (K 2 > 0) are shown in figures 6 and 7. In figure 6, hermaphroditism is the ESS. It can be seen from the plot, that the hermaphrodite population is robust against invasion by large numbers of females. Simulations started with low population densities evolved to the trivial equilibrium (not shown).

These observations are in agreement with mean-field kinetics. In figure 7, coexistence of females and hermaphrodites is evolutionary stable. The paths shown in figure 7 converge on an equilibrium at approximately (0.77, 0.12), this corresponds to the stable node, (h 5 , f 5 ), in the mean-field analysis. This is difficult to see from figure 7 due to slow dynamics close to the equilibrium.

Simulations started with low population densities evolved to the trivial equilibrium (not shown).

Breaking the condition K 2 > S (equation ( 9)) with K 3 > K 1 does not result in extinction following female invasion, as predicted by mean-field analysis. Instead simulations display coexistence with increasingly cyclic behaviour as K 2 is decreased and/or K 3 increased. This result is not surprising given the observations for K 2 = 0 shown earlier (figure 5).

The frequency of females in natural gynodioecious populations is often used as an indicator of the genetic mechanism responsible for male sterility in a particular species [START_REF] Bailey | A field guide to models of sex-ratio evolution in gynodioecious species[END_REF]. For cases of nuclear male sterility where K 2 > S, spatial simulations arrive at equilibria which are fairly static in (h, f )-space. Thus, equilibrium frequencies of females in the population are easily calculated. Figure 8 plots the equilibrium female frequencies for spatial simulation and mean-field analysis as a function of

K 3 /K 1 , with K 2 = K 3 .
For each set of results K 1 was fixed while K 2 and K 3 were incremented. By testing several values of K 1 , the plots show that the relative magnitudes of K 1 , K 2 and K 3 are dominant in determining the behaviour of the system and not their absolute values. Though birth rates do need to be of sufficient magnitude to avoid extinction. Reproductive rates used for the simulations in figures 5-8 are larger than those use for meanfield dynamics (figures 1-4). This because the strictly local interactions of the simulation make the population more vulnerable to extinction at low birth rates [START_REF] Preece | Sustainability of dioecious and hermaphrodite populations on a lattice[END_REF].

Mean-field analysis predicted that nuclear male sterility could evolve if 8 shows that the condition for nuclear male sterility to evolve in a spatial simulation with (

K 3 > K 1 . Figure
K 2 = K 3 ) is approximately K 3 > 3K 1 .
This infers that nuclear male sterility is less likely to evolve in species with local pollination and seed setting, as is the case in the simulation, than in than in a well mixed/widely dispersing population. Results from lattice simulations are summarised in table 2.

Figure 8 also shows that female frequency is lower in the spatial simulation than predicted by mean-field analysis. Thus, it may be expected that gynodioecious plant species with local pollination and seed setting should have lower female frequencies in the wild compared to species with long range interactions. The difference between mean-field and simulation is likely to be due to local variation in the concentration of females. As offspring are placed on nearest neighbour sites females with easy access to pollen can-

Male Sterility

Parameters Evolutionary Outcome

Nuclear 3K 1 > K 3 & K 3 = K 2 Hermaphroditism is ESS Nuclear K 3 > 3K 1 & K 3 = K 2 Gynodioecy Cytoplasmic 1.5 > K 3 K 1 & K 2 = 0 Hermaphroditism is ESS Cytoplasmic 8 > K 3 K 1 > 1.5 & K 2 = 0 Gynodioecy Cytoplasmic K 3 K 1 > 8 & K 2 = 0
Female invasion leads to extinction not compensate for females being starved of pollen in regions of high female concentration. This would result in regions with high female concentration dying off due to a shortage of pollen.

In many hermaphrodite plant populations individuals are self-incompatible [START_REF] Igic | The distribution of plant mating systems: Study bias against obligately outcrossing species[END_REF], which is the focus of the model presented above. In other species some self-fertilisation is observed in addition to out-crossing [START_REF] Kohn | Sex-ratio, seed production, biomass allocation, and the cost of male function in cucurbita-foetidissima (cucurbitaceae)[END_REF][START_REF] Ashman | Determinants of sex allocation in a gynodioecious wild strawberry: implications for the evolution of dioecy and sexual dimorphism[END_REF]. Though selfing would alter the details of the model presented in this article, we would still expect the main conclusions of our investigation to stand. Namely, females require a greater fertility advantage to become established when interactions are local, and the equilibrium frequency of females is lower when interactions are local. This may account for observations of [START_REF] Ashman | Determinants of sex allocation in a gynodioecious wild strawberry: implications for the evolution of dioecy and sexual dimorphism[END_REF] that observed female frequencies are generally lower than those predicted by the equation of [START_REF] Lewis | Male sterility in natural populations of hermaphrodite plants[END_REF] in a gynodioecious species of wild strawberry.

As pointed out in the introduction, male sterility has also been observed to be controlled by more complex genetic interactions than studied here. For example male sterility due to recessive genes [START_REF] Lewis | Male sterility in natural populations of hermaphrodite plants[END_REF] and the evolution of restorer genes which counter male sterility [START_REF] Schnable | The molecular basis of cytoplasmic male sterility and fertility restoration[END_REF].

Spatial simulations similar to the one presented here could easily be applied to these systems.

Conclusions

The population dynamics for a stochastic spatial simulation and its meanfield approximation have been presented for a gynodioecious population where male sterility is conferred by a dominant nuclear allele. It was demonstrated that a female needs a much greater fertility advantage to become established in the stochastic spatial simulation when compared with the mean-field approximation. This result suggests that gynodioecy is less likely to evolve in plants with local pollination and seed setting, and may also contribute to the observation that a relatively small percentage of plant species exhibit gynodioecy [START_REF] Richards | Plant breeding systems[END_REF] and dioecy [START_REF] Dellaporta | Sex determination in flowering plants[END_REF].

The equilibrium frequency of females in the population was found to be lower in the stochastic spatial simulation than in the mean-field approximation. This suggests that a lower frequency of females should be expected for gynodioecious populations with local pollination and seed setting.

Acknowledgements

We would like to thank the Royal Society, the EPSRC and the University of Nottingham for financial support. For helpful discussions we would like to thank J. P. Garrahan and A. Davison.

A. Derivation of mean-field dynamics

According to the reproductive process described in section 2, a hermaphrodite produces offspring via its female sexual function at a mean rate,

r = κ 1 4 4 j=0 [jP (j) P (x = |j)] , (13)
where j is the number of nearest neighbour sites occupied by hermaphrodites, κ 1 is the birth rate, P (j) is the probability a hermaphrodite has j nearest neighbours, and P (x = |j) is the probability that the site chosen for the offspring to reside on, x, is empty given j nearest neighbour sites are occupied by hermaphrodites. Applying Bayes' theorem to equation( 13),

r = κ 1 4 4 j=0 [jP (x = ) P (j|x = )] . (14)
As P (x = ) is not a function of j it can be taken out of sum in equation( 14)

and placed as a prefactor. In the meanfield,

P (x = ) = (1 -h -f ), (15)
where h and f are the densities of hermaphrodites and females respectively.

In the meanfield, the remaining sum from equation ( 14) is,

Figure 1 :Figure 2 :Figure 3 :Figure 4 :Figure 5 :Figure 6 :Figure 7 :

 1234567 Figure 1: Typical mean-field dynamics for the case where hermaphroditism is the ESS. Filled circles mark stable nodes, the empty circle marks a saddle point. The arrows indicate the direction of evolution. Plot shown for K 1 = 5, K 2 = 3 and K 3 = 3.

Figure 8 :

 8 Figure 8: Comparison of equilibrium sex-ratio for mean-field (line) and lattice simulation. Simulations shown for K 2 = K 3 , and K 1 = 10 ( ) and K 1 = 20 (×). Points shown were averaged over the last 10 generations of the simulation.

Table 2 :

 2 Evolutionary outcomes for lattice simulations, where K 1 , K 2 and K 3 are the scaled reproductive rates for the following processes: H + H → H, H + F → H, and H + F → F respectively (H represents a hermaphrodite and F a female). A value of K 1 > 6 is required for a sustainable hermaphrodite population.

The factor of 3 arises because P (j|x = ) = 0 for j = 4. Substituting equations( 15) and ( 16) into equation( 13) gives,

Equation( 17) is the meanfield reproductive rate per hermaphrodite via its female sexual function. A similar derivation can be used to calculate the meanfield rates that hermaphrodites and females are produced by female individuals.

B. Stability analysis -nuclear male sterility

The equilibrium points for nuclear male sterility are given in equations

(3) -( 8). The eigenvalues for the trivial equilibrium (0, 0) are,

with corresponding eigenvectors,

The eigenvalues for (h 2 , 0) are,

with corresponding eigenvectors,

The eigenvalues for (h 3 , 0) are,

with corresponding eigenvectors,

The eigenvalues of (h 4 , 0) are,

where

The corresponding eigenvectors are,

The eigenvalues of (h 5 , 0) are,

where the corresponding eigenvectors are,

C. Stability analysis -cytoplasmic male sterility

(h 2 ,0) and (h 3 ,0) are still equilibria for the case of cytoplasmic male sterility (K 1 > 0, K 2 = 0 and K 3 > 0). Though, (h 4 ,f 4 ) and (h 5 ,f 5 ) are not. For the special case K 3 = K 1 , a curve of equilibrium points exists given by equation ( 12). The eigenvalues from the linear stability analysis of points on this curve are given by,

where the hermaphrodite density h refers to the hermaphrodite density on the equilibrium curve. The corresponding eigenvectors are,

again the hermaphrodite density h refers to the hermaphrodite density on the equilibrium curve.


Dr Amke Caliebe

Arne Jochens

Michael Krawczak

Uwe Rösler
email: roesler@math.uni-kiel.de

A Markov chain description of the

Keywords: microsatellite, evolution, population genetics, asymptotic behaviour

come   L'archive ouverte pluridisciplinaire

Introduction

Microsatellites are successive iterations of a given short DNA sequence motif (usually 2-6 nucleotides long) that is repeated 5-100 times [START_REF] Tautz | Notes on the definition and nomenclature of tandemly repetive DNA sequences[END_REF][START_REF] Chambers | Microsatellites: consensus and controversy[END_REF]. The number of iterations (the "repeat number") serves to identify a given microsatellite allele. Microsatellites are abundant in many species and have very high mutation rates (up to 10 -2 per generation, [START_REF] Li | Microsatellites: genomic distribution, putative functions and mutational mechanisms: a review[END_REF]. Owing to their high degree of variability, microsatellites are frequently used as markers in population genetics [START_REF] Goldstein | The use of microsatellite variation to infer population structure and demographic history in a natural model system[END_REF][START_REF] Kashi | Simple sequence repeats as advantageous mutators in evolution[END_REF], DNA fingerprinting [START_REF] Cassidy | DNA testing in animal forensics[END_REF][START_REF] Bindu | Allele frequency distribution based on 17 STR markers in three major Dravidian linguistic populations of Andhra Pradesh, India[END_REF], whole genome mapping [START_REF] Weissenbach | A second-generation linkage map of the human genome[END_REF] and genetic epidemiology [START_REF] Thibodeau | Microsatellite instability in cancer of the proximal colon[END_REF][START_REF] Ashley | Trinucleotide repeat expansion and human disease[END_REF].

The stepwise mutation model (SMM) was first introduced by [START_REF] Ohta | A model of mutation appropriate to estimate the number of electrophoretically detectable alleles in a finite population[END_REF] to describe the behaviour of electrophoretically detectable alleles in a population. Since then, the SMM has been widely used for modelling microsatellite mutation and evolution [START_REF] Tishkoff | Global patterns of linkage disequilibrium at the CD4 locus and modern human origins[END_REF][START_REF] Zhivotovsky | Features of evolution and expansion of modern humans, inferred from genomewide microsatellite markers[END_REF][START_REF] De Iorio | Stepwise mutation likelihood computation by sequential importance sampling in subdivided population models[END_REF][START_REF] Vardo | Clonal diversity of a lizard malaria parasite, Plasmodium mexicanum, in its vertebrate host, the western fence lizard: role of variation in transmission intensity over time and space[END_REF]. The SMM assumes that, in one generation, the repeat number can only increase or decrease by at most one, usually with equal probability. More refined models have been proposed that include mutations of greater length, mutation rates that depend upon repeat number, or the additional introduction of point mutations [START_REF] Di Rienzo | Mutational processes of simple-sequence repeat loci in human populations[END_REF][START_REF] Garza | Microsatellite allele frequencies in humans and chimpanzees, with implications for constraints on allele size[END_REF][START_REF] Feldman | Microsatellite genetic distances with range constraints: analytic description and problems of estimation[END_REF][START_REF] Zhivotovsky | Biased mutations and microsatellite variation[END_REF][START_REF] Kruglyak | Equilibrium distributions of microsatellite repeat length resulting from a balance between slippage events and point mutations[END_REF][START_REF] Durrett | A new stochastic model of microsatellite evolution[END_REF][START_REF] Falush | Size-dependent mutability and microsatellite constraints[END_REF][START_REF] Calabrese | Dynamics of microsatellite divergence and proportional slippage/point mutation models[END_REF]; for an overview, see [START_REF] Watkins | Microsatellite evolution: Markov transition functions for a suite of models[END_REF] or [START_REF] Calabrese | Models of Microsatellite Evolution[END_REF]. As yet, however, it has remained controversial to what extent these models approximate the reality [START_REF] Chambers | Microsatellites: consensus and controversy[END_REF][START_REF] Whittaker | Likelihood-based estimation of microsatellite mutation rates[END_REF][START_REF] Sainudiin | Microsatellite mutation models: insights from a comparison of humans and chimpanzees[END_REF][START_REF] Cornuet | Inference on microsatellite mutation processes in the invasive mite, Varroa destructor, using reversible jump Markov chain Monte Carlo[END_REF].

In the following, we will consider the classical SMM. In 1975, Moran discovered that the distribution of the absolute frequencies n i (t) of alleles (as identified by their repeat number i) at time t does not converge, but has bounded variance.

He subsequently conjectured that the distribution "remains in a bunch" and characterised its behaviour as "wandering", without being more specific as to the existence of a limiting distribution [START_REF] Moran | Wandering distributions and the electrophoretic profile[END_REF]. To investigate convergence, Moran considered quantities C k (t) := N -2 i n i (t)n i+k (t), where N is the population size. For k = 0, this is the "effective number of neutral alleles in the population" of [START_REF] Ohta | A model of mutation appropriate to estimate the number of electrophoretically detectable alleles in a finite population[END_REF]. Moran was able to show that "unlike most problems in population genetics that have been discussed in the past, we do not obtain a limiting distribution or convergence in probability [of [START_REF] Moran | Wandering distributions and the electrophoretic profile[END_REF]. Shortly after Moran's publication, Kingman investigated the normalised Markov chain of the SMM, given by the repeat number difference to the allele of the N-th (or any other) individual in each generation [START_REF] Kingman | Coherent random walks arising in some genetical models[END_REF]. Using characteristic functions, he could prove exponentially fast convergence in distribution for a generalised model. He also obtained results about the limiting distribution of samples from a population when the population size tends to infinity conditioned that a certain relationship between time and population size holds.

C k (t)]."
Here, we will give a detailed analysis of the behaviour of the allele process under the SMM, where our focus will be upon the resulting Markov chain. Markov processes have been applied before to the characterisation of microsatellite mutation models by [START_REF] Watkins | Microsatellite evolution: Markov transition functions for a suite of models[END_REF]. In contrast to [START_REF] Kingman | Coherent random walks arising in some genetical models[END_REF], who used the analytic tool of characteristic functions, we will apply the stochastic method of recurrence of Markov chains. In Section 2 we will introduce the stepwise mutation model which is the basis for all subsequent results. In Section 3 the allele process X is investigated. We will make use of the fact that every population which does not die out, such as under a Wright-Fisher model, contains a genealogical lifeline that does not die out. Adding independent mutations to the genealogy generates an inherent random walk, and thereby results for the marginal distribution of X. In the second subsection, we will show that X is an irreducible, aperiodic and null recurrent Markov chain. The behaviour of X represents the global aspect of the SMM.

The normalised allele process V is analysed in Section 4 characterising the local view of the SMM. Again, marginal results such as moments and exponential moments will be given in the first subsection. Then, it will be proven that V is a positive recurrent Markov chain with exponentially fast convergence to the invariant distribution. A central result is provided in the third subsection where it will be shown that the marginal invariant distribution is unimodal. Finally, some simulation results for this distribution are given.

Wright-Fisher Model with Stepwise Mutations

The microsatellite allele process under neutral evolution will be studied using a

Wright-Fisher model with stepwise mutation. Let (Ω, A, P) be the underlying probability space, and let N ∈ N := {1, 2, 3, ...} be the constant population size. A microsatellite allele will be represented by the number of iterations of the sequence motif, the repeat number. Alleles are normalised such that allele 0 corresponds to a particular basic state m ∈ N, e.g. the most commonly observed repeat number. For simplicity in the classical SMM, which we apply here, there are no length restrictions on the allele size and even negative repeat numbers are theoretically possible. Thus, the set of possible alleles equals Z, and an allele z ∈ Z then has repeat number m + z.

a) Genealogy

The genealogy is assumed to be given by a Wright-Fisher model. For ease of notation and terminology, we will consider only haploid individuals. However, our results can easily be transferred to diploid individuals by regarding each of their two alleles separately.

Z n : Ω → {0, 1, -1} N is assumed to satisfy P(Z n (i) = 0) = 1 -μ, P(Z n (i) = 1) = P(Z n (i) = -1) = μ/2.
As usual, we assume an Independence Property for the genealogical and mutational processes, namely that the whole family

Z n (i), Y n (i), n ∈ N, i ∈ {1, ..., N} is independent. (1) c) Allele process
Let X n (i) denote the allele of the ith individual in the nth generation. For all n ∈ N 0 , X n : Ω → Z N can be written as

X n (ω)(i) := X n-1 (ω)(Y n (ω)(i)) + Z n (ω)(i) with X 0 ≡ 0 . (2)
X := (X n ) n∈N 0 is called the allele process of the Wright-Fisher SMM. The distribution of the initial states X 0 is arbitrary and does not influence the asymptotic behaviour. For the sake of simplicity, we assume X 0 ≡ 0 which means that all alleles have the same repeat number m at time 0.

d) Fundamental properties of the allele process X

Let A n := σ(Y 1 , ..., Y n , Z 1 , ..., Z n ) be the σ algebra generated by Y 1 , ..., Y n , Z 1 , ..., Z n .
Then the following property follows directly from the definition of X.

Proposition 1

(i) X n is A n -measurable for all n ∈ N. (ii) For all n ∈ N the family (X n (i)) i∈{1,...,N } is exchangeable. (Exchangeability Property)
3 Global Behaviour: the Allele Process X

Marginal Properties of X

To investigate the marginal distribution of the allele process X, we will use an immanent random walk. This is generated by the "lifeline" of the genealogy, i.e. the line of descent that never dies out. J n is the index, in generation n, of the (unique) member of the lifeline.

Proposition 2

(i) There exists an almost surely unique

J : Ω → {1, ..., N} N such that Y n (J n ) = J n-1 for all n ∈ N, and J n is σ (Y k , k ∈ N, k > n) measurable. Furthermore,
for n ∈ N, X n (J n ) has the same distribution as X n (1).

(ii) (X n (J n )) n∈N 0 is a random walk. For k ∈ Z, the transition probabilities are

P(X n (J n ) = k|X n-1 (J n-1 ) = k) = 1-μ and P(X n (J n ) = k-1|X n-1 (J n-1 ) = k) = P(X n (J n ) = k + 1|X n-1 (J n-1 ) = k) = μ/2.
Proof. (ii) follows from the definition of X once the existence of J has been established. Let τ n be the first generation (after n) in which all individuals have a common ancestor in generation n, i.e.

τ n := inf{k > n : ∀i, j ∈ {1, ..., N} Y n+1 • ... • Y k (i) = Y n+1 • ... • Y k (j)} . τ n is σ (Y k , k ∈ N, k > n) measurable and almost surely finite.
Then, for n ∈ N 0 , define on τ n < ∞

J n := Y n+1 • Y n+2 • ... • Y τn-1 • Y τn (1) .
J n is almost surely well defined. For

τ n = τ n-1 the equality Y n (J n ) = J n-1 is clear. For τ n > τ n-1 define Z := Y τ n-1 +1 • ... • Y τn (1). Then Y n (J n ) = Y n • Y n+1 • ... • Y τ n-1 (Z) = J n-1 .
Hence J n satisfies the required properties. 2

The first and second moment of the marginal distribution of X and a recurrence equation follow immediately from this proposition and from Prop. 1(i).

A limit result for the first absolute moment can be derived by applying the central limit theorem to the random walk (

X n (J n )) n∈N 0 . For all n ∈ N, z ∈ Z define ρ n (z) := P(X n (i) = z) .
Note that, owing to the exchangeability property of Prop. 1(ii), ρ n (z) is independent of the choice of i ∈ {1, ..., N}.

Lemma 3 For any

i ∈ {1, ..., N}, n ∈ N (i) E (X n (i)) = 0 (ii) Var (X n (i)) = μ n (iii) ρ n (z) = (1 -μ)ρ n-1 (z) + μ/2 (ρ n-1 (z -1) + ρ n-1 (z + 1)) (iv) lim m→∞ E |X m (i)| = ∞ Note that lim n→∞ Var (X n (i)) = ∞. Lemma 4 For any i, j ∈ {1, ..., N}, i = j, n ∈ N Cov (X n (i), X n (j)) = μ n + (N -1) n -N n N n-1 .
A proof of Lemma 4 is given in the appendix.

Characterisation of X as a Markov chain

The following theorem shows that, in our new representation as a Markov chain, the allele process X is null recurrent (see [START_REF] Breiman | Probability. SIAM[END_REF], p. 140, for the definition of null recurrent). Therefore, no asymptotic distribution exists.

In the following, we will write 0 N for (0, ..., 0) ∈ Z N . For the definition of A n , see Prop. 1.

Theorem 5

(i) The allele process X is an irreducible, aperiodic Markov chain on Z N with respect to (A n ) n∈N 0 .

(ii) The allele process X is null recurrent.

Proof. (i) follows directly from the definition of X. For the proof of the recurrence, it suffices to verify recurrence for state 0 N ∈ Z N because of irreducibility.

Remember that X 0 ≡ 0 N . We will prove the criterion , 1967, p.23, Theorem 4). One possibility for process X to get from state 0 N at time 0 to state 0

∞ n=1 P (X n = 0 N ) = ∞ (Chung
N at time 2n + 1, is that X 2n (1) = 0, Y 2n+1 (i) = 1
and Z 2n+1 (i) = 0 for all i ∈ {1, ..., N}. Therefore,

P (X 2n+1 = 0 N ) ≥ P (X 2n (1) = 0) 1 N N (1 -μ) N .
Now choose J according to Prop. 2(i). Then

∞ n=1 P (X n = 0 N ) ≥ ∞ n=1 P (X 2n (J 2n ) = 0) 1 N N (1 -μ) N = ∞
since the random walk of Prop. 2(ii) is known to be recurrent.

Let τ := inf{n ∈ N : X n = 0 N }. For null recurrence, it remains to be shown Definition 6 The process V := (V n ) n∈N 0 , defined by

that E (τ ) = ∞. This follows from τ ≥ inf{n ∈ N : X n (J n ) = 0}
V n : Ω → Z N -1 with V n (i) := X n (i) -X n (N) ,
is called the normalised allele process.

Marginal properties of V

In this subsection several marginal properties of V are derived. The proofs are given in the appendix.

The first and second moments of the marginal distribution of V can be calculated directly from the corresponding moments of X (see appendix). Because of the exchangeability property, the distribution of V n (i) is symmetric around zero.

Lemma 7 For any i, j ∈ {1, ..., N -1}, i = j, n ∈ N

(i) E (V n (i)) = 0, (ii) Var (V n (i)) = 2μN 1 -1 -1 N n , (iii) Cov (V n (i), V n (j)) = 1 2 Var (V n (i)).
Note that, in contrast to the behaviour of process X (see Lemma 3),

lim n→∞ Var (V n (i)) = 2μN is finite.
We now derive a recursion for the marginal distribution of V . Note that, because of the exchangeability property of Proposition 1(ii), the distribution of V n (i) is independent of i for i ≤ N -1. Thus, define

η n (z) := P(V n (i) = z) for all n ∈ N 0 , z ∈ Z . (3)
For k ∈ Z let r(k)

:= P (Z n (1) -Z n (2) = k) . (4)
Obviously r does not depend on n, r(0

) = 1-2μ+ 3 2 μ 2 , r(1) = r(-1) = μ-μ 2 , r(2) = r(-2) = 1
4 μ 2 and r(k) = 0 for any other k.

Lemma 8 For any

n ∈ N, z ∈ Z η n (z) = N -1 N 2 k=-2 r(k) η n-1 (z -k) + 1 N r(z) .
The next lemma provides a recursion for the higher moments and allows determination of the exponential moments of V n (i). (i) All moments of V n (i) are finite and emerge from the following recursion:

E (V n (i)) m = 0 for odd m. E (V n (i)) m = 1 N 2μ + μ 2 2 m-1 -2 + 1 - 1 N m k=0 k even m k 2μ + μ 2 2 m-k-1 -2 E V n-1 (i) k for even m.
(ii) All exponential moments of V n (i) are finite and are given by

E exp (λV n (i)) = 1 N n-1 k=0 1 - 1 N k c(λ) k+1 + 1 - 1 N n c(λ) n = c(λ) N 1 -1 -1 N n c(λ) n 1 -1 -1 N c(λ) + 1 - 1 N n c(λ) n .
The following corollary is straightforward and reveals the behaviour of the moments of V n (i) for n → ∞.

Corollary 10 Let i ∈ {1, ..., N -1}, m ∈ N, λ > 0.

(i)

lim n→∞ E (V n (i)) m < ∞. (ii) lim n→∞ E exp (λV n (i)) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ < ∞ if 1 -1 N c(λ) < 1 = ∞ if 1 -1 N c(λ) ≥ 1 .
Thus, process V fulfills the Doeblin condition and sections (ii) and (iii) follow (see, e.g. Doob (1953), pp. 192 ff., case b). 2

Unimodality of the asymptotic marginal distribution of V

Theorem 11 implies that the distribution η n of V n (i) (see Eq. ( 3)) converges in distribution as n → ∞. Let η = lim n→∞ η n . We will now show that η is a unimodal discrete distribution, which is one of our main novel results.

Following [START_REF] Keilson | Some results for discrete unimodality[END_REF], we call a distribution p on Z unimodal,

if at least one M ∈ Z exists such that p(z) ≥ p(z -1) for all z ≤ M p(z + 1) ≤ p(z) for all z ≥ M .
For proving the unimodality of η we need the following preparatory lemma, the proof of which can be found in the appendix.

Lemma 12 Let R + denote the set of strictly positive real numbers and R + 0 the set of positive real numbers including zero. If

M := {ν : Z → R + 0 |∃ n ∈ N; a 1 , . . . , a n ∈ R + ; b 1 , . . . , b n ∈ N 0 : ν = n i=1 a i • 1 {-b i ,-b i +1,...,b i } }, then M is closed under convolution.
With this, we can show that η is unimodal. The critical assumption of the following theorem, namely that μ ≤ 0.8, can safely be assumed for microsatellites.

Theorem 13

(i) If r : Z → R is defined as in Eq. ( 4), then η 1 = r and for all n ∈ N\{1}

η n = 1 N r * n-2 i=0 N -1 N i r i + N -1 N n-1 r n , ( 5
)
where * denotes the convolution of two functions and r i the ith convolution of r.

(ii) If μ ≤ 0.8, then η is unimodal and symmetric around zero.

Proof. Recalling that X 0 ≡ 0 N , it follows that, for all z ∈ Z,

η 1 (z) = P(V 1 (1) = z) = P(X 1 (1)-X 1 (N) = z) = P(Z 1 (1)-Z 1 (N) = z) = r(z),
according to the definition of r, see Eq. ( 4). Since by definition k∈Z r(k)η n-1 (z-k) = (r * η n-1 )(z), we can reformulate the recursive equation in Lemma 8 as follows:

η n = N -1 N r * η n-1 + 1 N r. ( 6
)
We will now prove Eq. ( 5) by induction. For n = 2, Eq. ( 5) follows from Eq. ( 6). Now, let n ∈ N\{1, 2}. Assuming that Eq. ( 5) holds for n -1, we have

η n = N -1 N r * η n-1 + 1 N r = 1 N r + N -1 N r * 1 N r * n-3 i=0 N -1 N i r i + N -1 N n-2 r n-1 = 1 N r * 1 + n-3 i=0 N -1 N i+1 r i+1 + N -1 N n-1 r n = 1 N r * n-2 i=0 N -1 N i r i + N -1 N n-1 r n .
To prove section (ii), we will first show that η n is unimodal for all n ∈ N.

Since 0 < μ ≤ 0.8, the following inequalities hold:

1 -2μ + 3 2 μ 2 > μ -μ 2 ≥ 1 4 μ 2 .
Thus, in the notation of Lemma 12, r

∈ M with n = 3, b 1 = 0, b 2 = 1, b 3 = 2, a 1 = 1-2μ + 3 2 μ 2 -(μ -μ 2 ), a 2 = -1 4 μ 2 +(μ -μ 2 ), a 3 = 1 4 μ 2 . Now, considering
Eq. ( 5), Lemma 12 implies that η n ∈ M for all n ∈ N, and all elements of M are clearly unimodal. Unimodality of η follows from the fact that the limit of a convergent sequence of unimodal discrete distributions is itself unimodal (see "Statement 4" in [START_REF] Keilson | Some results for discrete unimodality[END_REF]). 2

Simulation of the asymptotic marginal distribution η

Lemma 8 can be used to simulate the marginal distribution of V n (i). From

Theorem 11, we know that V n (i) converges in distribution as n → ∞. Figs.

1 and 2 show the behaviour in time of the distribution of V n (i), assuming μ = 0.01 and either N = 100 or N = 1000, respectively. For N = 100, the distribution of V n (i) is close to the stationary distribution at n = 100, and the domain is mainly concentrated in the interval [-7, 7]. For N = 1000, convergence is slower and the domain is larger. The distribution of V n (i) is close to the stationary distribution at n = 1000, and the domain is mainly concentrated in [-13, 13].

Discussion

We have shown that the allele process of the stepwise mutation model is characterised by two different types of behaviour. The expectation of the absolute value of the repeat number of a given individual converges to infinity. This signifies the global behaviour, where no convergence occurs. However, when the allelic state of an individual is chosen as a reference point for the other individuals of the population, then a limiting invariant distribution of the resulting allele difference process emerges. This is the local behaviour of the allele process, which implies that the alleles stay "clumped together" during convergence to infinity. These results confirm Moran's notion of the term "wandering distributions" [START_REF] Moran | Wandering distributions and the electrophoretic profile[END_REF].

The convergence of the allelic differences is exponentially fast, as was already noted by Kingman [START_REF] Kingman | Coherent random walks arising in some genetical models[END_REF]. This is reassuring because it means that estimates or test statistics obtained from the allele differences not only approach a limiting distribution, but do so very quickly. As we showed, the resulting limiting marginal distribution is unimodal.

It should be noted that the SMM is a very simple model of microsatellite mutation. In some cases, it would be reasonable to assume not only mutations that change the repeat number by one unit but to allow a wider range of mutations [START_REF] Huang | Mutation patterns at dinucleotide microsatellite loci in humans[END_REF]. Kingman also considered generalised forms of mutations [START_REF] Kingman | Coherent random walks arising in some genetical models[END_REF]. As long as the individual mutation events Z n (i) remain independent, which is biologically plausible, central Theorem 11 of this paper will hold true. If the random walk corresponding to the mutation process Z is null recurrent, Theorem 5 will apply. Another limitation of the SMM is the unboundedness of the state space whereas, in reality, negative repeat numbers cannot occur. Also, very large repeat numbers can result in physically unstable microsatellites and stop the evolutionary process at certain thresholds. One way to account for these limitations would be to restrict the state space of the allele process X by reflecting boundaries. The result would be a Markov chain with finite state space, and convergence to an invariant distribution would follow even for the non-normalised process X. However, differences between the normalised and non-normalised behaviour of the allele process remain possible, for instance, in the form of different convergence rates or different shapes of the invariant distribution. Because of the Markov structure and the assumed one-unit-up-or-down mutations, the process would only "realise" the existence of boundaries when it would be very close to them.

Most of the time, the process would stay away from the boundaries and behave according to the stationary distribution of the normalised process V , as if no boundaries would exist.

Therefore, from Lemma 14,

Cov (X n (i), X n (j)) = μ N N -1 N n + n(1 -N -1 N ) -1 (1 -N -1 N ) 2 = μ n + (N -1) n -N n N n-1 . 2
Proof of Lemma 7

(i) follows directly from the definition of V . Variance and covariance can be derived using Lemmata 3 and 4:

Var (V n (i)) = Var (X n (i)) + Var (X n (N)) -2Cov (X n (i), X n (N)) = 2μn -2μ n + (N -1) n -N n N n-1 Cov (V n (i), V n (j)) = Var (X n (N)) -Cov (X n (N), X n (j)) -Cov (X n (i), X n (N)) + Cov (X n (i), X n (j)) = Var (X n (N)) -Cov (X n (N), X n (j)) 2

Proof of Lemma 8

Using recursion (2), it follows that

η n (z) = P (X n (1) -X n (N) = z) = P (X n (1) -X n (2) = z) = P (X n-1 (Y n (1)) + Z n (1) -(X n-1 (Y n (2)) + Z n (2)) = z) = P (Y n (1) = Y n (2)) • P (X n-1 (Y n (1)) + Z n (1) -X n-1 (Y n (2)) -Z n (2) = z |Y n (1) = Y n (2)) + P (Y n (1) = Y n (2)) • P (Z n (1) -Z n (2) = z |Y n (1) = Y n (2)) = N -1 N k∈Z P (Z n (1) -Z n (2) = k) • P X n-1 (Y n (1)) + Z n (1) -X n-1 (Y n (2)) -Z n (2) = z Y n (1) = Y n (2), Z n (1) -Z n (2) = k + 1 N P (Z n (1) -Z n (2) = z) = N -1 N k∈Z r(k)P (X n-1 (1) -X n-1 (2) = z -k) + 1 N r(z) = N -1 N k∈Z r(k)η n-1 (z -k) + 1 N r(z) . 2

Proof of Lemma 9

(i): From Eq. ( 2), we obtain

E (V n (i)) m = P (Y n (i) = Y n (N)) E (Z n (i) -Z n (N)) m + N k,l=1 k =l P (Y n (i) = k, Y n (N) = l) • E (X n-1 (k) + Z n (i) -X n-1 (l) -Z n (N)) m = 1 N E ((Z n (i) -Z n (N)) m ) + 1 - 1 N E ((V n-1 (i) + Z n (i) -Z n (N)) m ) = 1 N E ((Z n (i) -Z n (N)) m ) + 1 - 1 N m k=0 m k E (Z n (i) -Z n (N)) m-k E V n-1 (i) k . Using E ((Z n (i) -Z n (N)) m ) = 0 for m odd and E ((Z n (i) -Z n (N)) m ) = 2μ + μ 2 (2 m-1
-2) for m even, section (i) follows by induction and Lemma 7.

(ii): Treating the exponential moments in the same way yields Since all the characteristic functions in the last expression are symmetrical around zero, it follows that ν * ν ∈ M. 2

d n := E (exp (λV n (i))) = 1 N E (exp (λ (Z n (i) -Z n (N)))) + 1 - 1 N E (exp (λV n-1 )) E (exp (λ (Z n (i) -Z n (N)))) = 1 N c(λ) + 1 - 1 N c(λ) d n-1 . 2 Proof of

Let

  Y n (i) be the direct ancestor of the ith individual in the nth generation, i ∈ {1, ..., N}, n ∈ N. Clearly, Y n (i) is an individual of the (n -1)th generation. According to the Wright-Fisher model, Y n : Ω → {1, ..., N} N satisfies P(Y n (i) = j) = 1 N for all i, j ∈ {1, ..., N} and the (Y n (i)) i∈{1,...,N },n∈N are independent. b) Mutation process Let Z n (i) be the mutational event preceding inheritance, from Y n (i), of the allele of the ith individual in the nth generation. We only consider mutation events that either increase or decrease the repeat number by 1, or leave the repeat number unchanged, i.e. Z n (i) ∈ {-1, 0, 1}. Let μ ∈ (0, 1) be the mutation rate, i.e. the probability of a change in repeat number per generation and per individual. Then,

4

  and from the fact that the random walk of Prop. 2(ii) is null recurrent. 2 Local Behaviour: the Normalised Allele Process V Since no asymptotic distribution exists for the allele process X, we will now consider the normalised allele process V , corresponding to the differences between the repeat numbers of each allele and the allele of the N-th individual in each generation. Note that because of the exchangeability, any other individual may take the place of the N-th individual.

Lemma 9

 9 For λ > 0 and i ≤ N -1, define c(λ) := E (exp (λ (Z n (i) -Z n (N)))). Then c(λ) = r(0) + (exp(λ) + exp(-λ)) r(1) + (exp(2λ) + exp(-2λ)) r(2). Let i ∈ {1, ..., N -1}, n, m ∈ N, λ > 0.

  Lemma 12Let * denote the convolution of two functions and define f b := 1 {-b,-b+1,...,b} for b ∈ N 0 . First note that, for all b ≤ b ∈ N 0 and for all z ∈ Z, the following equation holds:(f b * f b )(z) = i∈Z f b (i) • f b (zi) b +i,...,b +i} (z).Dropping argument z and decomposing the sum on the right-hand side,f b * f b = f b + -b +i,...,b +i} + 1 {-b -i,...,b -i} ) = f b + b i=1 (f b +i + f b -i ). Now let ν, ν ∈ M, which can be written as ν = n i=1 a i •f b i and ν = n i=1 a i •f b i .Let m(ij) := max{b i , b j }. Taking into account the distributivity and linearity of the convolution, the above implies thatν * ν = j • (f b i * f b j ) j • f m(ij)-k .

Figure 1 :Figure 2 :

 12 Figure 1: Convergence of the marginal distribution of the normalised allele process V .

Acknowledgments

This work was partly funded by the German Ministry of Science and Education (BMBF) through an NGFN SMP-GEM grant to Michael Krawczak (01GS0426).

Characterisation of V as a Markov chain

Like the original allele process X, the normalised process V is a Markov chain.

Contrary to X, however, V can be shown to be positive recurrent (see below;

for the definition of positive recurrent see [START_REF] Breiman | Probability. SIAM[END_REF], p. 140). Therefore, there is an invariant distribution that characterises the asymptotic behaviour of V , and V can even be shown to converge to this distribution exponentially fast. It should be pointed out that, whereas our Markov chain characterisation of the normalised allele process V is new, the convergence result was already obtained by Kingman, using characteristic functions [START_REF] Kingman | Coherent random walks arising in some genetical models[END_REF].

Theorem 11

(ii) V is positive recurrent.

(iii) V converges exponentially fast to the unique invariant distribution.

Proof. Using Eq. (2), section (i) follows from the fact that

For the proof of (ii) and (iii), write 0 N -1 for (0, ..., 0) ∈ Z N -1 and note that, for every z ∈ Z N -1 ,

Regarding the simplicity of the SMM, our results are only a first step towards a better understanding of the real-life situation, and investigations of how the allele process behaves under more realistic models incorporating, for example, variable mutation rates or migration, are warranted.

Appendix

In order to calculate the respective covariances of Lemma 4 from Eq. ( 2), we need the following little lemma that can be proven by induction.

Lemma 14 Let a, b ∈ R, b = 1. Then, for the real valued sequence (x n ) n∈N 0 defined by x n = (n -1)a + bx n-1 and x 0 = 0,

Proof of Lemma 4

Note that the exchangeability property implies that P(X n (i) = y, X n (j) = z) is independent of i, j ∈ {1, ..., N} as long as i = j. We can calculate a recursion for the covariances using Eq. (2), Prop. 1(i), Lemma 3(i), (ii) and the independence property (1):

y,z∈Z yz P(X n-1 (k) = y, X n-1 (l) = z) = 1 N Var(X n-1 (1)) + N -1 N Cov(X n-1 (i), X n-1 (j)) = 1 N μ(n -1) + N -1 N Cov(X n-1 (i), X n-1 (j)) .


Erik Hollnagel

From protection to resilience: Changing views on how to achieve safety

Effective safety management requires the ability to learn from the past and to anticipate the future. Yet what we can learn from the past (i.e., accident investigation) and what we can imagine for the future (i.e., risk assessment) depends critically on how we think about it, i.e., the models and methods we have at our disposal. Accident investigations have long been dominated by a search for causes, either as root causes or human errors. Risk assessment has similarly been dominated by static representations such as event and fault trees. In both cases the commonly used models and methods have reached their limits because the reality of our self-created socio-technical environments has become too complex. The alternative is to understand how the variability of human actions is a resource rather than a threat and to define safety as a system's resilience, its ability to adapt and adjust, rather than as the absence of adverse outcomes.

Introduction

Even in the best of all possible worlds, the future is not completely predictable. Events are bound to occur for which we are not prepared, some with positive and some with negative outcomes. Although there are very few situations where things go wrong compared to the very many where things work out fine and where the outcomes are as intended -or at least acceptable under the circumstances -the positive cases tend on the whole to go unnoticed. When the outcome of a task or an activity is acceptable, there is little motivation to look for why that was so; it is simply taken for granted -and even considered normal -that things go right. Conversely, when something goes wrong a relentless hunt for the cause(s) begins, in order to ensure that such an event never happens again. Unless we are willing to treat adversity with Panglossian optimism we must, of course, find some way to reduce the uncertainty, especially with regard to things that can go wrong. One venue for that is to design processes, systems, and organisations such that hazards are eliminated, or to ensure that the probability of adverse events is reduced to an acceptable level. (A hazard is here defined as an event that can lead to a known loss.) In order to do so, it is necessary that the system can be described in detail and that events develop in a predictable manner. However, since a growing number of sociotechnical systems are intractable, it is in practice impossible to achieve an acceptable level of safety by precautionary measures alone, i.e., by eliminating hazards, by preventing unexpected event, or by protecting against unwanted outcomes. Safety by design (analytical safety) must therefore be complemented by safety by management (operational safety). As resilience engineering makes clear, safety is something a system does, rather than something a system has [START_REF] Hollnagel | Resilience engineering: Concepts and precepts[END_REF].

Safety by design

The meaning of "safety by design" in this paper is that all possible, or practicable, precautions needed to ensure an acceptable level of safety are taken ahead of time. This could be either when the system is conceived of or designed, when detailed plans for operation are made, or when it is made ready for operation. The engineering parts of the system, namely the technology or hardware, are designed in detail and must perforce be configured, operated, and maintained according to meticulously prepared instructions. In such cases "safety by design" is not only an option but a requirement. Even so, there are limitations to what this can accomplish because technology always is a means to an end rather than an end in itself:

"Clearly not all accidents can be prevented by design. There are some consequences of technology that cannot reasonably be predicted at the design stage, particularly in new technologies, using new materials and scientific principles. However, once these have led to accidents, there is a clear responsibility for designers to prevent them in future designs. " Hale, Kirwan, Kjellen (2007, p. 310) Yet even if the technology is reliable, and even if it functions in a highly predictable manner, there are other problems. For instance, have the safety requirements for the detailed design been adequately defined? Are the selected concepts proven from a safety point of view? Or will will it be possible to implement regulatory, corporate and customer safety requirements within acceptable cost limits? Although these questions address issues that have to do with the use of the system in the wider operational context, they are legitimate for "safety by design" as well. The non-technical parts of the system, namely the people or the liveware, are far more difficult to design. Indeed, an organisation can never be explicitly designed in the same way that a piece of machinery can, one fundamental reason being that the "components" by their nature are variable and flexible, regardless of whether they are considered individually or collectively. Humans and social systems simply do not function like machines, despite courageous attempts to make them do so through training, interaction design, and automation. For the non-technical parts the alternative is to focus on the operation of the system, and on how it is possible to keep the variability within acceptable limits, i.e., by managing safety.

Safety by management

Safety management has during the last decades of the 20th century become a major issue in itself, with a number of commercial solutions for safety management systems (SMS) on offer. It may therefore be sensible to look beyond the use of SMS as a glib phrase and try to understand what it entails. The crucial point is here the definition of safety, i.e., the definition of what the SMS is supposed to manage. Safety management is, from a practical point of view, a kind of control; indeed, the meaning of 'manage' is 'to exert control over, to direct.' Safety management can therefore be interpreted as the control of the organisational functions and practices that together produce safety. The purposes of an SMS is to ensure that the organisation's "safety processes" develop in the indended direction and that they are not disrupted or hindered by internal or external events and conditions. This reformulation leads to three subsidiary questions. The first is what the organisational "safety processes" really are, i.e., what it is that "produces" safety. The second is how these processes can be controlled in practice, i.e., how their "speed" and "direction" can be changed. And the third is how the outcome or the result can be measured, i.e., what the proper indicators of safety are. As far as the latter is concerned, it is common practice to associate safety with the "freedom of unacceptable risk," or as ICAO (2006) defines it:

"Safety is the state in which the risk of harm to persons or of property damage is reduced to, and maintained at or below, an acceptable level through a continuing process of hazard identification and risk management." By defining safety in this way, a safe state is defined by the absence of something and the measures or indicators are therefore the, hopefully, diminishing numbers of negative events. The established practice is to try to reach this state through a tripartite approach: to eliminate hazards (by design), to prevent initiating events that may lead to adverse outcomes (by constraining operations), and to protect against adverse outcomes (by introducing barriers). This approach to safety implies a distinction between normal and abnormal operations. The normal operations ensure that the (safety) process goes in the right direction and that it produces what it was supposed to produce or achieve. The abnormal, or off-normal, operations disrupt or disturb the normal operations or otherwise render them ineffective. The main motivation for safety management, as it is commonly practiced, has therefore been to prevent such disruptions or disturbances from taking place. In that sense safety management has been driven by what has happened in the past, hence been mainly reactive.

Theory W

One way of characterising the established approaches to safety and safety management is to propose two idealised positions, expressed as two different theories called Theory W and Theory Z, respectively. According to Theory W, socio-technical systems are safe and efficient because:

• The systems are well-designed and scrupulously maintained.

• The procedures that are provided to operate the systems are complete and correct.

• The operators of the systems, the people at the so-called sharp end, behave as they are expected to, and as they have been trained to.

• Designers can foresee every contingency and therefore provide the systems with appropriate response capabilities. Theory W describes a world of well-designed, well-tested, and well-behaved systems. It is a characteristic of such systems that there is a high degree of reliability of equipment; that workers and managers are vigilant in their testing, observations, use of procedures, and operations; that staff is well trained; that management is enlightened, and that good operating procedures are in place. The threat to normal performance comes from different types of failures or malfunctions, such as active failures and latent conditions, equipment faults, and human error. Safety can therefore be achieved by constraining performance variability in various ways (Figure 1).

It follows from the assumptions behind Theory W that performance variability of any kind should be avoided. Technological systems are designed to perform one -or a few functions -very efficiently and with limited variability. Functions that depend on humans, whether as individuals, as groups, or as organisations, are more versatile but also less uniform and the variability is seen as a threat or a disturbance that may end in performance failures. Everything possible is therefore done to prevent this from happening. The change from normal to abnormal operations can happen either abruptly, as when something breaks, or gradually in the form of a drift or slow misalignment [START_REF] Cook | Going solid": a model of system dynamics and consequences for patient safety[END_REF]. The solution according to Theory W is therefore to constrain performance, for instance, by means of barriers, interlocks, rules, procedures, standardization, interaction design, norms, etc., or even by replacing humans by technology, as in the use of automation. This is supplemented by more wide-ranging activities such as risk management (comprising, e.g., hazard identification, risk assessment, risk mitigation, and risk communication), hazard and incident reporting, safety investigations, safety analyses and safety studies, and safety performance monitoring.

Safety for underspecified systems

Theory W corresponds to systems where the principles of functioning are clear and where there is good predictability. This was not an unreasonable assumption at the time when the currently established methods were developed, roughly between 1965 and 1985, but it is much less reasonable today. There are two main reasons for that, first that all systems of interest are more or less intractable, and second that performance variability is inevitable. In order for a system to be tractable, four conditions must be fulfilled: (1) that the principles of functioning are known, (2) that a description does not contain too many details, (3) that a description can be made relatively quickly, and (4) that the system does not change while the description is being made. The last condition is the most important one, and is in a way a synthesis of the three preceding. Theory W clearly requires that systems are tractable. Many of the present day systems of major interest for industrial safety are unfortunately intractable rather than tractable. This means that the principles of functioning are only partly known, that the description is elaborate and contains many details, that it takes a long time to make, and that the system therefore change while the description is made. In consequence of that it will be impossible to provide a complete description or specification of the system. Intractable systems are underspecified in the sense that details may be missing or unavailable (e.g., [START_REF] Clarke | Safety culture: Underspecified and overrated?[END_REF]. But if a system is underspecified it is clearly not possible to provide precise procedures or instructions. The people working in the system, be it at the sharp end or at the blunt end, must therefore be able to apply the available prescriptions and procedures to conditions and situations that differ from what was assumed. In other words, it is necessary that people are able to vary or adapt what they do to ensure that the system functions as required and achieves its operational goals. Performance variability -whether it is called improvisation, adaptation, efficiencythoroughness trade-off, sacrificing decisions, or creativity -is therefore necessary, hence an asset rather than a threat.

The inevitability of performance variability

While machines and technological artefacts are designed, built, and maintained so that they can produce a near constant performance -at least until they fail and must be replaced -the same is not the case for humans and for organisations. There are many reasons why human performance never can be constant or machine-like. One is that physiological functions (muscles, nerve cells, sensory organs, etc.) are subject to fatigue, saturation, and accommodation, and that they regularly require a period of rest or reconstitution. Many basic psychological functions, such as attention or vigilance, are also limited with regard to how long they can be maintained at a constant level. A second reason is that humans seem to have an innate tendency to vary what they do, often to avoid monotony or to find an easier way to accomplish a task. It is this ingenuity and creativity that is at the heart of adaptability and of the ability to overcome constraints and underspecification. A third reason is socially induced variability, for instance in the sense that others will have expectations -or even informal normsabout how much and how little effort is acceptable in a given situation. A different source of variability is the organisational culture, specifically the safety culture. A fourth reason is that performance depends on the state of the organisation and environment, which may vary in terms of demands to work, resources available, etc. A fifth and final reason is the variability due to ambient working conditions, for instance temperature, climate, humidity, noise, etc. Performance variability is for these reasons inevitable, whether on the level of the individual or the organisation. This is fortunately not a disadvantage since performance variability is necessary to overcome the underspecification of large socio-technical systems. But it is a problem if safety models and methods fail to recognise that.

Theory Z

Since most systems of interest today are intractable, it is impossible to provide a complete description of them or to specify what an operator should do even for many normally occurring situations. This is recognised by Theory Z, according to which socio-technical systems are safe and efficient for the following reasons:

• Humans learn to overcome the inevitable shortcomings, such as design flaws and functional glitches.

• Humans can adjust their performance to meet the actual demands of a situation.

• Humans can interpret procedures and apply them to suit actual conditions.

• Humans can detect when somethings fails or goes wrong, and can in many cases correct for it as well.

In Theory Z, performance variability is both normal and necessary and is the source of positive and negative outcomes -successes and failures -alike. Failures can consequently not be prevented by eliminating performance variability, i.e., safety cannot be managed by constraints. The solution is instead to identify the situations where normal performance variability may combine to create unwanted effects and to monitor continuously how the system functions in order to intervene and "dampen" performance variability that threatens to get out of control, cf. Figure 2. Conversely, performance variability should be accentuated or amplified when it can lead to successful outcomes. Thus rather than looking for ways in which something can fail or malfunction, we should try to understand the characteristics of normal performance variability, and specifically how internal and external factors may affect the size and nature of the variability.

Safety by management: Resilience engineering

The difference between Theory W and Theory Z corresponds to the difference between a reactive and a proactive approach to safety management. If performance variability is both normal and necessary, safety must be achieved by managing performance variability rather than by constraining it. This is consistent with resilience engineering, which is based on the following principles [START_REF] Hollnagel | Resilience engineering: Concepts and precepts[END_REF].

• Performance conditions are always underspecified, as argued above, and individuals and organisations must therefore adjust their performance to the current conditions. Because resources and time are finite, such adjustments will inevitably be approximate. • For tractable systems, most adverse events can be attributed to a breakdown or malfunctioning of components and normal system functions. For intractable systems, most adverse events cannot. They are instead best understood as the result of unexpected combinations of normal performance variability or as the converse of the adaptations necessary to cope with real-world complexity. • Effective safety management cannot be based on hindsight, nor rely on error tabulation and the calculation of failure probabilities. Safety management cannot be only reactive but must also be proactive. Resilience Engineering looks for ways to enhance the ability of organisations to create processes that are robust yet flexible, to monitor and revise risk models, and to use resources proactively in the face of disruptions or ongoing production and economic pressures. A resilient system is defined by its ability to adjust its functioning prior to, during, or following changes and disturbances so that it can go on working even after a major mishap or in the presence of continuous stress. A resilient system accepts a constant sense of unease and remains sensitive to the possibility of failure [START_REF] Hollnagel | Remaining sensitive to the possibility of failure[END_REF]. The quality of resilience can be defined more precisely by pointing to four essential qualities or abilities that a system or an organisation must have, cf. Figure 3.

• A resilient system must be able to respond to regular and irregular threats in a robust, yet flexible, manner. It is not enough to have a ready-made set of responses at hand, since actual situations often do not match the expected situations -the only possible exceptions being routine normal operation. The organisation must be able to apply the prepared response such that it matches the current conditions both in terms of needs and in terms of resources. In terms of the three types of threats proposed by [START_REF] Westrum | A typology of resilience situations[END_REF], this is the ability to deal with regular threats. The responses enables the organisation to cope with the actual.

• A resilient system must be able flexibly to monitor what is going on, including its own performance. The flexibility means that the basis for monitoring must be assessed from time to time, to avoid being trapped by routine and habits. The monitoring enables the system to cope with that which is, or could become, critical in the near term.

• A resilient system must be able to anticipate disruptions, pressures, and their consequences. This means the ability to look beyond the current situation and the near future, and to consider what may happened in the medium-to long-term. In terms of the three types of threats proposed by Westrum (op. cit.), this is the ability to deal with the irregular threats, possibly even the unexampled events. The anticipation enables the system to cope with the potential.

• Finally, a resilient system must be able to learn from experience. This sounds rather straightforward, but a concrete solution requires consideration of which data to learn from, when to learn, and how the learning should show itself in the organisation -as changes to procedures, changes to roles and functions, or changes to the organisation itself. The learning enables the organisation to cope with the factual.

Conclusions

Effective safety management requires the ability to learn from the past and to anticipate the future. Yet what we can learn from the past (i.e., accident investigation) and what we can imagine for the future (i.e., risk assessment) depends critically on how we think about it and on the models and methods we have at our disposal. Accident investigations have long been dominated by a search for causes, either as root causes or human errors. Risk assessment has similarly been dominated by representations of linear combinations of a small number of events, such as event and fault trees. In both cases the underlying assumption is that systems -and events -are tractable. It makes sense that models and method would be just about adequate for the typical type of problems at the time they were developed. Indeed, there would be little reason to develop a method that was more complex or more powerful than required, not least because it would be difficult to imagine what that should comprise. New models and methods are developed because existing models and methods sooner or later encounter problems for which they are inefficient or inadequate. This, in turn, happens because the socio-technical systems where accidents happen continue to develop and to become more complex and more tightly coupled. The inevitable result is that any method after a while becomes underpowered because the nature of the problems change, although it may have been perfectly adequate for the problems it was developed for in the first place. The risks that dominate in present day systems have a different aetiology than the risks that dominated even one or two decades ago. This has two important ramifications. The first is that it is more difficult to understand present day risks -at least until an accident has happened. It is harder to understand the "mechanisms," because risks can arise from non-linear interactions among normal performance variability as well as from consequences of failures and malfunctions. And because of that it is also more difficult to think of ways to reduce or eliminate the risks, hence to manage safety. The second ramification is that many of the established risk assessment and accident investigation methods are inadequate for tightly coupled, intractable systems. This dilemma was made clear when [START_REF] Perrow | Normal accidents: Living with high risk technologies[END_REF] proposed that accidents could be seen as normal, in contrast to risk assessment and accident investigation methods that naturally focus on that which is abnormal or dysfunctional. The lesson to be learnt from that is that we must continue to evaluate critically the methods that are at our disposal. The fact that a method has worked in the past is no guarantee that it will also work in the future. Indeed, by the time a method has become adopted as a standard it is almost certain to be outdated. The ways that socio-technical systems develop means that risks can emerge in different ways, and that existing methods therefore sooner or later will need to be complemented with more powerful approaches. What these will be, no one can say for certain.

Figure 1 .

 1 Figure 1. Safety according to Theory W

Figure 2 :

 2 Figure 2: Safety according to Theory Z

Figure 3 :

 3 Figure 3: Principles of resilience engineering