{'fraction_non_alphanumeric': 0.03907253907959953, 'fraction_numerical': 0.007173418811867207, 'mean_word_length': 4.427498467198038, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Prediction of invasion from the early stage of an epidemic

Predictability of undesired events is a question of great interest in many scientific disciplines including seismology, economy, and epidemiology. Here, we focus on the predictability of invasion of a broad class of epidemics caused by diseases that lead to permanent immunity of infected hosts after recovery or death. We approach the problem from the perspective of the science of complexity by proposing and testing several strategies for the estimation of important characteristics of epidemics, such as the probability of invasion. Our results suggest that parsimonious approximate methodologies may lead to the most reliable and robust predictions. The proposed methodologies are first applied to analysis of experimentally observed epidemics: invasion of the fungal plant pathogen \emph{Rhizoctonia solani} in replicated host microcosms. We then consider numerical experiments of the SIR (susceptible-infected-removed) model to investigate the performance of the proposed methods in further detail. The suggested framework can be used as a valuable tool for quick assessment of epidemic threat at the stage when epidemics only start developing. Moreover, our work amplifies the significance of the small-scale and finite-time microcosm realizations of epidemics revealing their predictive power.


Introduction
Predictability of catastrophic events such as earthquakes, epidemics, fracture or financial crashes [1][2][3] is a topic of increasing interdisciplinary interest. The predictability of these events is inextricably linked to the inherent complexity of the phenomena under consideration [1,2]. Here, we focus on epidemiology. Within this context, many studies have been devoted to prediction of the temporal incidence of epidemics (i.e. the evolution of the number of infected hosts in the course of an epidemic) [4][5][6][7][8]. Recently, an increasing number of papers have also considered the prediction of the spatio-temporal evolution of epidemics [9][10][11][12].
Both the temporal and the spatio-temporal incidence depend on complex factors related to the transmission of infection and the properties of the hosts. For instance, the hosts are not identical in susceptibility and transmissibility of infection due to difference in age, size, genotype and neighbourhood. [5,9,10,13]. The transmission of infection is stochastic, meaning that a healthy host is infected by contact with inoculum from an infected host with a certain probability only. Many epidemics, notably those involving transmission by invertebrate vectors, or by wind and rain for many plant pathogens, are subject to variability in weather. This environmental stochasticity can also influence the evolution of epidemics in such heterogeneous systems [14]. All these factors make prediction of disease incidence an extremely challenging and sometimes controversial task [4,15,16]. In Ref. [4], it was suggested that, although obtaining precise quantitative predictions for the incidence would be obviously desirable, qualitative predictions may be more valuable. This is very much along the lines of ideas from the science of complexity claiming that, despite the fact that giving accurate predictions for the detailed evolution of complex systems might be an illusory task, certain qualitative features of the evolution, such as occurrence or absence of a catastrophic event, could be more amenable for prediction [2,17,18].
Here, we address the question of predictability of epidemics using a methodological framework inspired by the science of complexity. The main aim is to estimate the probability that an emerging epidemic will 2 invade a significant fraction of the population in the future. This quantity can be viewed as a qualitative feature of the complete spatio-temporal evolution of epidemics.
We propose several methods for approaching the problem that offer different levels of precision. Our results suggest that the most precise methods do not necessarily lead to more reliable predictions. Instead, parsimony seems to be the key ingredient for prediction based on inherently limited observations. The framework presented below deals with epidemics caused by a broad class of pathogens leading to permanent immunity of infected hosts after recovery (or death). There are numerous examples of such diseases affecting populations of humans [19,20], animals [21] and plants [22]. The advantage in analysis of such epidemics is that they are characterised by a well-defined final state consisting of only hosts that were never infected and hosts that were infected and became immune. In particular, we focus on the estimation of the probability that an epidemic is invasive in the final state.
The proposed methods are first applied to prediction of invasion of a pathogen in an experimental model system in which the fungal plant pathogen, Rhizoctonia solani, spreads through a population of hosts represented by discrete nutrient sites. The properties of the sites, e.g. nutrient concentration, can be varied for different realisations of epidemics. Such a system is convenient for generation and observation of rapid, highly replicated and repeatable epidemics and it is used as a benchmark for description of our methodologies and analyses. The epidemic prediction analysis for experimental system is followed by a test of our methods in numerical experiments for epidemics spreading on networks of hosts arranged on a regular lattice. The advantage of investigating such epidemics is that their properties are known beforehand and this allows us to provide a precise analysis of the performance of prediction methods by comparing with the expected behaviour.

Methods
The methods follow four steps that are basic for any scientifically meaningful prediction of the behaviour of a complex system: (i) Observation of the initial evolution of the process over a certain period of time, t obs . (ii) Construction of a model for description of the observed behaviour. (iii) Fitting the model to the data obtained from observations. (iv) Extrapolation of the behaviour to the future by using the model with the fitted parameters. These steps are interconnected and we argue that they should be kept at a similar level of complexity in order to make their interplay as consistent as possible. In this paper, we investigate whether or not such consistency is important for obtaining reliable predictions by exploring several combinations of strategies for steps (i)-(iii) (see summary in Table 1).
For concreteness, we illustrate our prediction methods for the particular case of fungal colony invasion in microcosms comprising populations of nutrient sites (agar dots) [23]. In these experiments, the central agar dot in ensembles with the geometry shown in Fig. 1 was inoculated by the soil-borne fungal plant pathogen R. solani and the spread of the fungal colony is scored in discrete time steps (e.g. daily).

Methodology
Step (i): data Step (ii): model Step (iii): fitting Parameters Table 1. Overview of methods for prediction. Combination of possible data sets considered in step (i) and methods used for addressing steps (ii) and (iii) in the prediction process. The first column introduces a label for each set of methods, which are ordered according to their overall expected precision. The second column gives the format of the data obtained from observations in step (i). The smallest precision of the data corresponds to methods A and B in which only the incidence, C(t), is used. Methods C-E use a limited spatio-temporal knowledge of the evolution of the infection given by the shell-evolution function F (l, t). Method F uses the time of infection of each host, {t i }, which is typically unknown from observations but can be inferred in step (iii). In step (ii), RF and CT are abbreviations for the Reed-Frost and continuous-time model, respectively. In step (iii), we have used several methods for fitting: Minimum Distance (MD), Approximate Bayesian Computation (ABC), and data-augmented Markov Chain Monte Carlo (DA-MCMC). Column five lists the parameters involved in each step for prediction of the fungal invasion in the agar-dot experiment. Methods based on the RF dynamics are parametrised by the transmissibility, T , and a time scale τ exp . The CT dynamics used in methods E and F is parametrised by T , a characteristic time τ 0 , and a shape parameter for the time-dependence of the transmission of infection, k. The RF dynamics corresponds to the limit k → ∞ of the CT dynamics. The MD method for fitting consists in minimising the parameter d 2 (last column) which measures the difference between observations and numerical simulations. The ABC fitting procedure assumes that a simulated invasion fits well the observed epidemic if d 2 < ǫ, where ǫ is a free parameter. As shown in the last column, the definition of d 2 depends on the descriptors for observations used in step (i). For methods A and B, d 2 is defined in terms of the observed and simulated incidences, c obs (t) = C obs (t)/N and c sim (t) = C sim (t)/N , normalised to the number of hosts in the population, N . In methods C-E, d 2 is defined in terms of the shell-evolution function for observations and simulations. Sections I and II of SI Appendix give more details on the definition of models used in step (ii) and fitting methods used in step (iii).
In the following, we give a general description of the steps for prediction of epidemic invasion with particular assumptions suitable for the analysis of the fungal invasion experiment. The main details of all the methodologies are summarised in Table 1. The methodology C is mainly used for illustration of our concepts. Details of other explored methodologies are given in SI Appendix. The motivation for choosing methodology C is twofold: (i) it keeps all the steps for prediction of the behaviour at a similar level of complexity, as illustrated by analysis of the fungal colony invasion in the agar-dot experiment; (ii) it is a parsimonious methodology that leads to predictions that are, at least, as robust as (or, arguably, even more robust than) those based on more sophisticated approaches (Table 1).
Step (i). The information that can be extracted from observation of the time evolution of epidemics is usually limited. In many cases, the only available information is the incidence, C(t) (the number of infected hosts), at subsequent observations, and occasionally the spatial location of infected hosts is also known, e.g. for epidemics in populations of plants [8,22,24,25]. These limitations have a dramatic influence on subsequent steps in the prediction process and it is crucial to identify which quantities are sufficient for prediction of the catastrophic event (i.e. the probability of an invasive epidemic in our case). As shown in Table 1, we consider three types of observations. The first consists of discrete temporal observations of C(t) (Methods A and B in Table 1). The second possibility (Methods C-E) considers discrete spatio-temporal observations giving the evolution of infection at discrete times, t, in shells at a 'chemical distance' l from the initially inoculated host. As explained in Fig. 1(a), such observations can be properly described in terms of a shell-evolution function F (l, t). As a third possibility (Method F), we use a method for data augmentation in step (iii) that infers the unobserved time of infection for each host, {t i } [8,24,25].
Step (ii). We describe the evolution of the epidemic in terms of a spatial SIR (susceptible-infectedremoved) epidemiological model where the hosts can be either susceptible (S), infected (I) or removed (R) [19,21,22,26]. This is a prototype model for a wide class of epidemics where disease leads to permanent immunity of hosts after recovery or death. In particular, this paradigm has been shown to be appropriate for description of fungal invasion [23,27,28]. In principle, a continuous-time dynamic model is necessary to provide a precise description of epidemic evolution characterised by stochasticity in times of infection and removal/recovery of hosts. Following this idea, it would be natural to use a model with continuous-time (CT) dynamics (described in SI Appendix ). The drawback of this approach is that it requires knowledge of the precise times of infection of hosts, {t i }, which are typically not available from discrete spatio-temporal observations in step (i). In order to match an appropriate model with the level of detail of observations, we consider the discrete-time dynamics model which reduces the SIR framework to the so-called Reed-Frost (RF) model [29]. This simplified description is not expected to capture the dynamical details of the evolution of the epidemic but reproduces well its final state [30,31]. This is a very important consequence of the fact that, no matter how complicated the evolution of the epidemic is, the final state of an epidemic with death of infected individuals or permanent acquired immunity after recovery depends only on the probability T , called transmissibility, that the infection has ever been passed between each pair of connected hosts (as shown in Fig. 1(b)). Although the transmissibility is expected to exhibit a certain degree of spatial heterogeneity in real epidemics, we make the minimal assumption that the trend of the epidemic can be well approximated by a RF process with a homogeneous effective transmissibility T .
Step (iii). The goal of this step is to estimate the values of the parameters of the model used in step (ii) that give a good description of the observations. Consider, for definiteness, methodology C in Table 1. Due to factors such as stochasticity and heterogeneity in transmission, a given observed spatiotemporal map for infection can occur for different values of the estimated transmissibility, T . However, some of these values for T are more likely to produce the observed spatio-temporal pattern than others. To account for this, we introduce the probability density function (p.d.f.), ρ( T ), which quantifies the probability that the observed spatio-temporal pattern is reproduced by a certain value of T . As shown schematically in Fig. 1(c), ρ( T ) is calculated by generating a large number of stochastic realizations 5 a Step (i): observations b Step (ii): model c Step (iii): fitting d Step (iv): estimated P inv (L) Steps for prediction used in methodology C (Table 1). (a-left) Discrete spatio-temporal observations of evolution (spatio-temporal map) of a hypothetical epidemic spreading from the central host in a population of hosts (circles) arranged on a triangular lattice with lattice spacing a and L hosts per side of the hexagonal boundaries of the system. The arrangement of nutrient sites in the fungal invasion experiment analysed below is of this type. Numbers inside the circles denote the times, t = 0, 1, . . ., of infection of hosts by the observation time t obs = 2. Empty circles correspond to healthy (susceptible) hosts by the same time. The hexagons with dotted lines indicate the shells of hosts at a given chemical distance, l, to the centre of the system. (a-right) The spatio-temporal evolution of the epidemic is described by the shell-evolution function F (l, t) giving the relative number of hosts in layer l infected by time t. For instance, F (2, 2) = 3/12 for the epidemic shown in the (a-left) panel. (b) In step (ii), the epidemic is described in terms of an SIR model. An infected host, I , remains infectious during the infectious period τ which, for simplicity, is taken as being constant over the whole population and set as a unit of time, τ = 1. After the infectious period τ , the host is removed, R . During the time τ , I can transmit the infection to a neighbouring susceptible host, S , with probability T (transmissibility). Alternatively, I can be removed without passing the infection to S with probability 1 − T . In the discrete time Reed-Frost (RF) dynamics used in our approach, the infection is passed instantaneously from I to S at t = τ . (c) In step (iii), the fitting procedure consists in finding the probability density function (p.d.f.) ρ( T ) for transmissibilities T (c-left) such that a RF process with T and shell-evolution function F sim (l, t) give a good description of the observed shell-evolution function, F sim (l, t) (c-right). For visualization clarity in c-right, an example of F sim (l, t) represented by the blue-grid surface does not fit well the observed shell-evolution function, F obs (l, t) (the shaded surface). The p.d.f. ρ( T ) is obtained by running many RF epidemics with random transmissibility and minimising the parameter d that quantifies the difference between the observed and the RF shell-evolution functions. of the RF epidemic with transmissibilities T sampled uniformly from the interval [0, 1] and comparing their shell-evolution functions, F sim (l, t) (see caption to Fig. 1 for definition), with the observed one, F obs (l, t). Ideally, the distribution ρ( T ) would correspond to the histogram of values for T producing shell-evolution functions F sim (l, t) identical to F obs (l, t) but obtaining an exact match is computationally very demanding. Moreover, reproducing F obs (l, t) by a RF model is in general impossible in realistic epidemics for which time is not discrete. Therefore, we use a minimum distance (MD) algorithm to calculate ρ( T ) approximately as the histogram of values of T minimising the quantity d 2 f defined in Table 1 that measures the distance between F sim (l, t) and F obs (l, t). In this way, the sampled values of T reproduce F obs (l, t) approximately rather than necessarily with distance d 2 f = 0. This approach is similar to the Approximate Bayesian Computation (ABC) method that determines ρ( T ) as the histogram of values of T for which d 2 f ≤ ǫ, where ǫ is a parameter used in the method [32]. Both ABC and MD algorithms give similar results despite the fact that MD does not require the use of an additional parameter ǫ. In addition to these approximate methods, we have fitted the spatio-temporal evolution of the CT model proposed in step (ii) for comparison by means of a more standard Bayesian procedure using Markov chain Monte Carlo (MCMC) method with data augmentation (method F in Table 1).
Step (iv). In the final stage of the prediction process, given ρ( T ), we evaluate the probability P inv (L) that the observed epidemic will ever invade a system of size L. The epidemic is defined as being invasive if the final cluster of removed hosts has reached at least one node on each of the six edges of the system. Otherwise, the epidemic is classified as being non-invasive. The conditional probability of invasion P inv ( T ; L) in a system of size L by an SIR process with a given transmissibility, T , can be calculated numerically by running many stochastic realisations of the epidemic and counting the fraction of invading events. As shown in Fig. 1(d), P inv ( T ; L) exhibits a sigmoidal dependence on T which indicates a non-invasive (invasive) regime of epidemics for relatively small (large) values of T . Once P inv ( T ; L) and ρ( T ) are known, the estimated probability of invasion can be calculated as follows: This formula defines the probability that the invasion occurs given our knowledge about the effective transmissibility encoded by ρ( T ) (see a simple graphical interpretation in terms of the shaded area in Fig. 1(d)). Importantly, Eq. (1) gives an extrapolation of the behaviour of the epidemic to its final state without necessarily providing a detailed description of the actual evolution leading to such a state.

Application to fungal invasion
In the fungal invasion experiments, the spatio-temporal maps of infected agar dots were scored daily over 21 days (see two typical patches of colonisation after 21 days in Fig. 2(a)). The transmissibility in this experiment corresponds to the probability of fungal colonisation between two adjacent agar dots and it was controlled by variable lattice spacing, a = 8, 10, 12, 14, 16, 18 mm. Clearly, the experimental setup is restricted both in space and time. Our aim is to use these limited observations to estimate the probability of invasive epidemics in larger systems and for longer times. The analysis is performed for each individual realisation of the experiment (6 replicates per value of a).
In order to make a proper comparison between the experimental observations with the RF model used in methods A-D, it is necessary to rescale the time step of the RF dynamics with dimensionless τ = 1 to τ exp measured in days. The value of τ exp is not known and it is treated at the same level as the transmissibility. More explicitly, we deal with a bi-variate probability density function, ρ 2 ( T , τ exp ), which can be determined for each epidemic with a simple extension of the methods explained in Sec. 2 (step (iii)) for obtaining ρ( T ). The estimated P inv is obtained from Eq. (1) by defining ρ( T ) as the marginal p.d.f.,  Fig. S1 in SI Appendix shows similar plots for six replicates of the system.
As explained in more detail in Sec. I of SI Appendix and summarised in Table 1, the continuous-time SIR model used in methods E and F involves three parameters: T ,τ 0 , and k. The fitting of the data results in a p.d.f. ρ 3 (T, τ 0 , k) from which we obtain ρ( T ) = ∞ 0 ρ 3 ( T , τ 0 , k)dτ 0 dk. The probability of invasion is then calculated from Eq. (1) in the same way as for methods A-D.

Uncertainty of the estimated transmissibility
The functions ρ( T ) obtained for the fungal invasion experiments typically exhibit a pronounced peak (see the results for one replicate in Fig. 3 and similar results for more replicates in Fig. S1 of SI Appendix ). The peaked shape of ρ( T ) suggests that T can be suitably described in terms of its mean value T and standard deviation σ T = ( T 2 − T 2 ) 1/2 . For each methodology, Fig. 4 shows the average over replicates of T and σ T as a function of the lattice spacing. These estimates correspond to observations of the evolution of infection during t obs = 21 days. All the methods give similar values for T which have a clear and expected tendency to decrease with increasing a. The uncertainty in T , quantified by σ T , exhibits greater variations between methodologies but it takes values that are smaller than T for all the methods and lattice spacings ( Fig. 4(b)). This means that T is a good measure of the typical value of the transmissibility. However, the value of T on its own does not necessarily provide a good approximation for P inv because the width of ρ( T ) can bring a significant contribution to the integral in Eq. (1). This is explicitly shown in Sec. 4.
Comparison of T , σ T , and ρ( T ) for different methods leads to the following conclusions: (1) Given a level of description (step (i)) and a model (step (ii)), the estimates of the transmissibility obtained with ABC and MD methods are, in general, in good agreement (cf. method A with method B, and method C with method D in Figs. 3 and 4).
(2) Given a level of description (step (i)) and an estimation method (step (iii)), the posteriors obtained  Table 1. using discrete-and continuous-time models are in reasonable agreement (cf. method D with method E in Figs. 3 and 4). The only difference is a trend for ρ( T ) corresponding to CT dynamics to have a "heavy tail" for large values of T (see the replicate 4 in Fig. S1, SI Appendix ). Large values of T are correlated with large values of the time scale τ 0 (i.e., slower processes with high T ) and small values of the shape parameter k, that are ruled out by the RF model. This effect becomes more important for larger values of the lattice spacing, as indicated by the large values of T and σ T corresponding to method E (asterisks in Fig. 4).
(3) The estimates from augmented-data MCMC (methodology F) are in general different from those obtained by other methods. Moreover, the p.d.f.'s ρ( T ) obtained with MCMC show no systematic trend with respect to the other methods. With respect to, e.g., the ρ( T ) obtained with MD, they can be located at slightly higher (replicates 5 and 6 in Fig. S1) or lower (replicates 2 and 3) values of T , or approximately at the same value (replicates 1 and 4). Moreover, the variation in the peak position of ρ( T ) between different replicates is larger than for the other methods. This suggests that the MCMC method is more sensitive to fine details of the evolution of the epidemic. A possible explanation is that augmented-data MCMC involves the inference of the unobserved colonisation times and thus is intrinsically individual-based, in contrast to shell-based (or mean-field) methods, which try to match the colonisation times in an approximate manner only.

Comparison of fitted models with experimental data
In order to assess the quality of the assumptions used for estimation, we compare the fitted models with the available experimental data. For methods A-D, we compare the incidence and shell-evolution function obtained numerically (RF dynamics) with values for T and τ exp sampled from ρ 2 ( T , τ exp ) estimated by means of the spatio-temporal maps at maximum observation time t obs = 21 days with the actual incidence and shell-evolution function for each epidemic. Similarly, the fits from methods E and F are compared with experimental observations by running numerical epidemics with parameters for the CT dynamics sampled from the p.d.f. ρ 3 ( T , τ 0 , k). We make a quantitative comparison based on squared distances d 2 c and d 2 f (cf. Table 1) between simulated epidemics and experimental fungal invasions. More explicitly, we define the root mean square (rms) distances, where ∆t = 21 days is the time interval used for calculations of d 2 c or d 2 f . The quantity l max is the maximum chemical distance to the centre of the system of agar dots. Its value decreases with the lattice spacing and ranges from l max = 2 for a = 16 mm to l max = 8 for a = 8 mm [23]. From the definition of d 2 c given in Table 1, it is easy to see that ∆c gives the typical deviation of the simulated incidence per unit host at a given time, c sim (t), from the observed incidence per host at the same time, c obs (t). Similarly, ∆F gives the typical deviation of the simulated shell-evolution function, F sim , evaluated at any spatio-temporal coordinates (l, t) from the observed value at the same coordinates. Fig. 5 shows the mean of the rms distances obtained by averaging over stochastic simulations and over replicates with given lattice spacing. The low values of the rms distances (∆c 0.2 and ∆F 0.3) indicate that the observed C(t) and F (l, t) are statistically well described by the fitted models. For any given method and lattice spacing, we obtained ∆c < ∆F , which is expected because reproducing the spatio-temporal evolution represented by F (l, t) is more demanding than capturing the temporal evolution of the colonisation given by c(t). Both ∆c and ∆F tend to be larger for a around 10-12 mm which, as shown below, corresponds to cases that are close to the invasion threshold (i.e. where P inv decreases from 1 to 0 on increasing a, as shown by the continuous line in Fig. 2(a)). Variability between replicates of epidemics with given T is larger around the invasion threshold, which is associated with a critical phase transition and characterised by large fluctuations [13, 21,23,[26][27][28]. As a consequence, the quality of fits is lower in the vicinity of the invasion threshold and this leads to larger values of ∆c and ∆F .
Methodology C gives a good balance between performance and number of parameters involved. Methods C-E based on an approximate spatio-temporal description of epidemics given by F (l, t) result in more accurate predictions than methods A (squares in Fig. 5) and B (diamonds) that neglect spatial features of invasion. Moreover, the approximate methods C-E also perform better than even methodology F (triangles in Fig. 5) despite the fact that the latter aims for a more precise spatio-temporal description. A more qualitative and visual comparison of estimated and observed C(t) reveals similar differences between all the methodologies (see details in Sec. III of SI Appendix ).

Two applications for prediction methods
As a first application of the proposed methods, we have studied the predictive power of the estimates of the probability of invasion and the incidence by calculating ρ( T ) from the early stages of the actual epidemic, i.e. for t obs < 21 days. In particular, based on the estimated ρ( T ) for t obs = 10 days, we have obtained estimates for the probability of invasion at time t = 21 days (squares and dashed line in Fig. 2(a)) and compared them with the probability P exp of invasion at t = 21 days obtained directly from the experimental data. The observed probability of invasion, P exp , is estimated by counting (for each a) the fraction of replicates in which the fungus has reached the six outer edges of the experimental system by 21 days. The mean of P inv averaged over replicates with the same value of a (symbols in Fig. 2(a)) gives a reasonable estimate for the observed mean of P exp (solid curve in Fig. 2(a)) after t = 21 days for most of the methods. Overall, the best predictions for P inv are obtained with methodology C (solid  Table 1, as marked in the legend. circles in Fig. 2(a)). As expected, P inv decreases with increasing a for all methodologies, illustrating the existence of the threshold for epidemics around a ≃ 12 mm [23]. Similarly, the experimentally observed incidence and shell-evolution function are statistically well captured by the numerical extrapolation for their simulated counterparts up to time t = 21 days obtained from observations over times t obs < 21 days. A visual illustration of the agreement between the observed and predicted incidence is given in Sec. III of Appendix SI. Fig. 6 shows the rms distances between the observed and predicted evolutions from day 11 to day 21. The time interval used to calculate ∆c and ∆F from Eqs. (2) and (3) is ∆t = 11 days. The relative trends of ∆c and ∆F between methods and lattice spacings are similar to those reported in Fig. 5 for the comparison between observations and numerical simulations with t obs . The main difference is that the values of the rms distances corresponding to predictions of the evolution of colonisation (Fig. 6) are systematically larger than those obtained by simply comparing observed evolutions with their respective fittings (Fig. 5). This is in agreement with the intuitive idea that predicting the a priori unknown evolution of a system is more challenging than reproducing a fitted evolution. As a second application of our methodology, we have calculated P inv at the end of the epidemic as a function of the lattice spacing in systems of size L = 51, i.e., larger than the experimental samples of sizes L = 2, ..., 8 which decrease with increasing lattice spacing (see the two populations for different value of a shown in Fig. 2(a)). Such predictions are based on estimates for the transmissibility obtained from observations up to day t obs . As expected, P inv decreases with increasing a, illustrating the existence of the threshold for epidemics around a ≃ 12 mm [23]. The results of applying each of the prediction methods are shown in Fig. 2(b) for the mean probability averaged over replicates for each value of a. All the methods except E give similar predictions for P inv . The large values of P inv predicted by method E are a consequence of the "heavy tail" of the p.d.f. ρ( T ) which gives a significant weigh to the high values of P inv for large T in Eq. (1). The dependence of P inv on a differs from the observed probability of invasion, P exp (continuous curve in Fig. 2(b)). The difference can be qualitatively understood by recalling that P inv gives an extrapolation both in space and time. Indeed, P inv ≥ P exp because some epidemics that are non-invasive after 21 days have a certain probability to invade a system of size L = 51 for  Table 1, as marked in the legend.
t > 21 days. In addition, both infectivity and susceptibility are expected to be subject to heterogeneity in the agar-dot system due to, e.g., inherent variability. Based on the results presented in Sec. 4 for numerical experiments with heterogeneity in transmission, we expect the estimated P inv to give an upper bound to the actual probability of invasion. This can also contribute to the difference between P inv and P exp for large values of a.

Numerical experiment
The quality of the predictions presented in the previous section is influenced by the quality of the observations in step (i), the suitability of the model chosen in step (ii) for description of the data, and the fitting procedure used in step (iii). In principle, the effect of these factors on predictions could be minimised by optimising the procedures used in each step for prediction. Stochasticity associated with the transmission of infection also influences the ability of making reliable predictions. In contrast to the previous factors, stochasticity is inherent to the nature of the system and its negative effect on predictions cannot be minimised without modifying the system. In this section, we present a sensitivity analysis of our methods by applying them to prediction of invasion for numerically-simulated epidemics where the main factor compromising predictability is the intrinsic stochasticity in transmission of infection. The advantage in this case with respect to more realistic situations is that both the transmissibility, T , and the probability of invasion, P inv (T ; L), are known and it is then possible to investigate the performance of the estimates for ρ( T ) and P inv (L) by comparing with the known quantities.
We first consider the simplest situation when the observed epidemics follow the RF dynamics with homogeneous transmission (i.e. T is the same for all pairs of nearest neighbours in the population). The idea is to run numerical experiments with known T , observe the evolution of the epidemic over an initial Hosts are placed on the nodes of a triangular lattice of size L = 51 (cf. Fig. 1(a)). The evolution of epidemics starting from the central host is observed over time t ≤ t obs = 7τ . The line marked by circles shows the dependence of the conditional probability of invasion, P inv (T ; L), on transmissibility obtained by the simulations in the system of size L = 51. The shaded region shows the levels of confidence in percentage of the p.d.f ρ( T ) around the most probable transmissibility, T ⋆ (continuous black line) corresponding to each value of P inv (L). The horizontal dashed line illustrates the case of estimations giving P inv (L) = 0.2. If the observed epidemic has a value of the transmissibility such as T l that is to the left of the curve for P inv (line with circles), the estimated P inv (L) overestimates P inv . In contrast, for values of the transmissibility that are to the right of the curve for P inv (e.g. T h ), the probability P inv (L) underestimates P inv .
interval of time, t obs , and then apply the methods described above to calculate P inv (L) assuming that T is unknown (as it occurs for real epidemics). For concreteness, we consider the arrangement shown in Fig.1(a) and use methodology C for steps (i)-(iii). RF epidemics are observed during t obs = 7τ with the aim of estimating the probability that they will invade a system of size L = 51. Note that over the time interval t ≤ t obs = 7 the epidemic at most invades a hexagon of size L = 15. Then, the behaviour of the epidemic is extrapolated both in space and time. We proceed by, first, calculating the p.d.f. ρ( T ) for the estimated transmissibilities T compatible with observation (spatio-temporal map). Fig. 1(c) shows an example of ρ( T ) obtained from the analysis of the evolution of an epidemic with T = 0.4. In general, the most probable estimate for the transmissibility, T * , corresponding to the maximum of ρ( T ) for a single epidemic, differs from T but not significantly. In many cases, T lies within the 68% confidence interval for ρ( T ) around its maximum (see a more detailed discussion in SI Appendix ). The distribution ρ( T ) allows the probability of invasion P inv (L) in the system of size L = 51 to be estimated using Eq. (1). We have applied this prediction method to many (∼ 10 4 ) spatio-temporal maps created with known transmissibility T spanning the interval [0, 1]. For each value of estimated P inv (L), the distribution ρ( T ) is represented by a horizontal slice of the shaded area in Fig. 7 (see e.g. the slice along the dashed blue line corresponding to P inv (L) = 0.2 with darker colour corresponding to higher probability relative to the maximum of ρ( T )). The black ridge in the shaded area corresponds to the most probable transmissibility, T ⋆ , for each value of P inv (L).
To test the quality of the predictions, the estimate of the probability of invasion P inv ( T ; L) is compared with the probability P inv (T ; L) that would be obtained if the exact value of T was known a priori (see line marked by circles in Fig. 7). Making such a comparison we can see that for epidemics with low transmissibility where invasion is possible but not highly probable, P inv (L) overestimates P inv (T ; L) for most of the possible values of T contributing to P inv (L) (the shaded area including the ridge region corresponding to the typical values of T is mainly above the line marked by circles in Fig. 7). This means that the estimations are biased upwards and the most likely is that the actual probability of invasion will be smaller than predicted. In other words, such predictions will typically give a safe bound for the probability of invasion. Obviously, there is a non-zero probability that the observed epidemic has a large value of the transmissibility (such as T h in Fig. 7). In this case, P inv (L) would underestimate the actual probability P inv . For more invasive epidemics (i.e. epidemics with P inv 0.5), the shaded area is mainly below the line marked by the circles in Fig. 7 meaning that the predicted probability of invasion P inv underestimates P inv for most of the possible values for T , including the most probable, T ⋆ . In these situations however, both P inv and P inv are large and the predictions allow for a reasonable assessment for invasion to be done. In Sec. VI of SI Appendix we show mathematically that differences between P inv (L) and P inv evaluated at the most probable transmissibility T ⋆ are mainly dictated by the curvature of P inv around T ⋆ and is intrinsically linked to the non-zero width of the p.d.f. ρ( T ). This general result implies that the biases in the probability of invasion at low and high transmissibility are independent of the fitting method (Table 1) because all methods lead to a p.d.f. ρ( T ) with non-zero width. The results presented above correspond to RF epidemics with homogeneous transmission. A similar approach has been used to deal with more realistic epidemics where the transmission of infection is heterogeneous due to variability in the infectivity and the susceptibility of hosts. As already mentioned in the previous section, the estimated P inv (L) for such epidemics usually gives a bound to the actual probability of invasion that is even safer than that obtained for cases with homogeneous transmissibility (see SI Appendix for more detail).

Discussion
The methodology introduced here focuses on the prediction of relatively simple but important features of epidemics. This is in contrast to much previous work dealing with the prediction of quantitative properties of epidemics such as the detailed spatio-temporal evolution of the incidence. The advantage in dealing with simple characteristics of epidemics is that they can be more easily predicted in terms of simplified description of the spatio-temporal evolution. Our results demonstrate that, under quite general assumptions, it is possible to give reliable prediction of the final state of an epidemic with permanent immunisation from the early stage of its evolution. Such a prediction is possible even for a single realisation of an epidemic and thus the framework is relevant to inherently unique real-world epidemics. In fact, our approach can be applied for prediction of epidemics in real systems characterised by a wide range of space and time scales (e.g. crops) based on micro-or meso-cosm experiments of finite size and over finite time.
The results obtained for experimental fungal invasion by using approximate methods (C-E in Table 1) are more robust than those based on supposedly more precise methodology F. This might be a consequence of the interplay between the high fitting precision for MCMC methods involving data augmentation with a poor description of the actual dynamics given by the continuous-time model fitted to the observations. A model capturing dynamical details at a level consistent with that offered by the fitting procedure might exhibit more predictive power than that presented here. This is an illustration of the importance of keeping all the steps involved in prediction at a similar level of complexity in order to give reliable predictions. For practical applications, the particular method to be used for obtaining the most reliable prediction depends on the problem in hand. The general rule that seems to emerge from our analysis is that a reasonable method should use as much informations as available from observations, and avoid inferring data that is not directly available unless this is strictly required by the problem. This is the case for methods C-E in our particular study of fungal colony invasion in the population of agar-dots. Methods A and B use less information than available from observations (i.e. they use C(t) instead of F (l, t)) while method F infers information that is not available from observations.
The proposed methods assume that epidemics can be approximately described by an effective transmissibility that is constant over time and homogeneous in space. However, their applicability goes beyond epidemics with constant transmissibility. In particular, we have shown that the method gives reliable predictions in the presence of spatial heterogeneity in the transmission of infection. We expect that the methododology can also be successfully applied to cases where the transmissibility changes over time but it remains within the bounds for the effective transmissibility estimated from the early stage.
We have characterised the final state of epidemics by the probability of invasion. This quantity is suitable for systems with well-defined boundaries such as, e.g. the population of agar-dots analysed above. In cases where hosts are placed on the nodes of more complex networks, the boundaries of the system are not necessarily well-defined [33] and it is more convenient to characterise the final state of epidemics in terms of the mean number of removed (i.e. ever infected) hosts, N R . Our approach also applies to such complex networks. The formula given by Eq. (1) provides an estimated size, N R (L), if P inv ( T ; L) is replaced by the function N R ( T ; L) giving the size of SIR epidemics with given transmissibility, T .
Our methods have been applied under the assumption that the network of contacts between hosts remains unchanged during the course of epidemics. Such an approximation has been widely used in the past [33] and it is reasonable for cases in which the rate of change of the configuration of contacts is much smaller than the removal rate of hosts (i.e. ∼ τ −1 in our notation). This condition is clearly satisfied for the fungal colony invasion of the population of agar-dots considered here and also for many other epidemics associated with pathogens spreading in, e.g., networks of plants [22], farms [9], or airports [11]. This paradigm is also applicable to the spread of many infections in human populations (for instance, measles or SARS that have recovery periods of the order of few days). In contrast, the dynamics of contacts between humans plays an important role for other infectious diseases such as syphilis with a recovery period of 100 days [34]. A possible strategy to make predictions in this kind of networks would involve inferring the mixing parameter for contacts (as defined, for instance, in [34]) based on observations (step (i)) in a similar way as we estimated the parameters of the SIR model in step (iii).
Another interesting task would be to extend the ideas presented here to deal with epidemics with persistence where immunity after recovery is not permanent (i.e. recovered hosts can be re-infected). In this case, the simplest model for description of observations is the SIS model (susceptible-infectedsusceptible) [35] and a possible quantity to be predicted would be the stationary prevalence of infection (i.e. the density of infected hosts in the stationary state reached after a transient [35]).
Due to stochasticity in transmission of infection, it is not possible to determine the parameters of a model describing an epidemic with absolute certainty even if the epidemic is observed during a long time t obs before attempting inference. Furthermore, if it were possible to determine the exact value of the parameters, it would still be impossible to make arbitrarily precise predictions of the evolution of the epidemic in the future or predict with absolute certainty if the epidemic is going to be invasive or not (instead, one has to deal with the probability of invasion). The uncertainty in the prediction of the evolution of epidemics grows monotonically with the look-ahead time (see, e.g. the forecast of the incidence in Sec. IV of SI Appendix ). There exists a prediction horizon beyond which the uncertainty of predictions of the evolution of the epidemic is too large for predictions to be useful. The location of the horizon is epidemic-dependent and also depends on how precise we want our predictions to be. In contrast, for a pure SIR epidemic, there is no prediction horizon for quantities such as P inv or N R that only depend on the transmissibility. In other words, different replicates of epidemics with given T will follow different evolutions that have a prediction horizon but will lead to the same P inv or N R [26,30,31]. More complex nonlinear dynamics for transmission associated with, for instance, a seasonal component in the transmission rate may lead to chaotic behaviour [36]. Predictability of catastrophic events in systems exhibiting chaotic behaviour is a non-trivial question that has been widely studied in the past [37] and still receives considerable attention at present [38]. Even in the absence of stochasticity, the prediction horizon in these systems is intrinsically limited due to the high sensitivity of chaotic processes to the initial conditions and the values of the parameters. In addition, the accuracy of predictions does not necessarily increase monotonically with the observation time, t obs , before prediction [39]. In such situation, it would be necessary to estimate the value of t obs leading to the most reliable prediction. Due to all these factors, the methods proposed in this paper may not work when applied for prediction of catastrophic events in nonlinear dynamical systems. However, the ideas presented here together with approaches proposed for prediction in nonlinear dynamical systems may help in devising strategies for prediction in stochastic non-linear systems.

I. CONTINUOUS-TIME EPIDEMIOLOGICAL MODEL FOR STEP (II)
In this section, we present the three-parameter continuous-time (CT) model used in methods E and F (Table 1, main text) to address step (ii) for prediction. In principle, this model is more realistic than the Reed-Frost (RF) model used in methods A-D and allows the existence of possible effects caused by the discrete-time character of the RF description to be explored. We use the fungal invasion experiment as a benchmark for the comparison between different models. In the CT, the spread (transmission) of the fungal colony between two neighbouring dots is treated as a time-inhomogeneous Poisson process [1].
The waiting time distribution f (t) for each transmission event can be modelled by a Weibull distribution multiplied by the transmissibility T : where τ 0 is the characteristic time scale of the process, and k is a shape parameter of the distribution.
Given the cumulative probability function P (t) = t u=0 f (u)du, the survival function S(t) (giving the probability that transmission did not occur by time t) obeys the following relation: The rate of the transmission process, φ(t), is a function of the time since colonisation of the donor dot and is given by the expression: In the limit k → ∞, the CT model reduces to the RF model, with the same value of T and infectious period τ = τ 0 . Indeed, for k → ∞, the survival probability given by Eq. (S.2) becomes a step function, which corresponds to a RF model, in which infection can be transmitted from an infected host to a susceptible neighbour with probability T only once the infectious period τ 0 has passed.

II. METHODS FOR FITTING MODELS TO DATA IN STEP (III)
The aim of this section is to give a detailed description of the different methods used for fitting models M to data D described in step (iii) of the methods for prediction proposed in the main text.
The ideas presented in the main text can be framed within a Bayesian approach which assumes that the parameters θ describing M are random variables. The aim of step (iii) is to evaluate π(θ|D), the probability density that M with parameters θ describes the data. According to Bayes' rule: where π(θ) is the prior distribution of the parameters, reflecting our initial belief in their values, and È(D|θ) is the likelihood (the probability of the data given the parameters).
Several challenges arise when using this Bayesian approach in the analysis of epidemic spread. The first difficulty is associated with inherent limitations in the observations. A complete spatio-temporal data set would contain the precise time of colonisation t j of each host j in the population, i.e. D = {t j }.
Unfortunately, it is often the case that observations do not provide such detailed information. In the particular experimental data set considered in the main text, the status of each dot is only recorded at discrete (1-day) time intervals. Hence, the actual dataset is D = {d j }, where d j is the day when dot j was first observed as colonised. We have explored two possible ways to deal with the lack of precise information. The option used in methods A-E (Table 1, main text) involves identifying descriptors for the evolution of the epidemic that are suitable for the prediction of the catastrophic event. The lowerdimensional descriptors used in this work are the cumulative incidence, C(t), and the shell-evolution function, F (l, t). Another option is to use data augmentation [2] that treats the unobserved colonisation times as parameters to be estimated. This is the procedure followed in method F.
Another source of difficulties is due to the fact that the analytical calculation of the posterior π(θ|D) is, in general, impossible [3]. Therefore, it is common to resort to numerical methods to sample from π(θ|D). In order to do this, we propose a new approximate method (denoted as MD), which is based on the calculation of the minimal distance between two datasets and does not require knowledge of the likelihood È(D|θ). In addition, we have used two known methods: • a method belonging to the class of Approximate Bayesian Computation (ABC) [4], which calculates an approximate posterior, and which shares several basic features with the MD method.
The ABC method was used with multiple purposes: (a) to test the new MD method (involving the minimisation of a given distance between observed and simulated data) against an already-known method (ABC, which involves a cutoff on the same distance).

A. Minimum distance (MD) method
Let D obs be the observed data, θ a set of candidate parameters, and D sim a simulated dataset generated using θ. Then, if D obs = D sim , the vector of parameters, θ, is drawn from π(θ|D obs ). In practice, obtaining an exact match between observed and simulated dataset is often computationally unfeasible, and one has to resort to an approximate match. To this end, we define a metric d 2 (D 1 , D 2 ) that measures the distance between two datasets D 1 and D 2 . The aim of the MD method is to calculate the distribution of parameters θ that minimise d 2 (D 1 , D 2 ) and gives an approximate posterior. The choice of the metric d 2 (D sim , D obs ) is problem-specific and, in general, not unique. We have tested two different metrics, corresponding to different descriptors of the data (cf. step (i) of the main text): 1. For the shell-based description, we used d 2 f = l, t (F sim (l, t) − F obs (l, t)) 2 , where l enumerates the shells, t is the discretised time (observation times in days), and F is the shell-evolution function.
2. Another option is to ignore any spatial information from the data (mean-field (MF) description) and consider only the total fraction of colonised sites at time t, c(t) = C(t)/N . In this case, we chose the distance function d 2 c = t (c sim (t) − c obs (t)) 2 , where t is again the discretised time.
The algorithm to implement the MD method considers two indexes, n, counting the parameter vectors resulting from the minimisation procedure, and r, counting the iterations. The values of n and r are in the range n, r ∈ N with maximum values n max and R, respectively, and proceeds as follows: For all the results presented in the paper and obtained with the MD algorithm, we set R = 5000.
For some fungal epidemics, we have checked that larger values of R do not lead to smaller values of

B. Approximate Bayesian Computation (ABC)
In common with the MD method, the ABC approximate Bayesian method we have used also relies on the definition of a metric, d 2 . If the distance between observed and simulated datasets is less than a given tolerance parameter, ǫ, i.e. d 2 (D obs , D sim ) ≤ ǫ, then θ is drawn from the approximate posterior π θ|d 2 (D obs , D sim ) ≤ ǫ . The accuracy of the approximation increases as ǫ → 0.
For our estimations, we used a Markov-chain Monte Carlo algorithm that can be summarised as follows: ABC.1 Set n = 0 and choose the initial value θ 0 of the parameter vector.
ABC.6 Set θ n+1 = θ ′ with probability p acc , or θ n+1 = θ n with probability 1 − p acc . ABC.7 Set n = n + 1 and return to ABC.2 until the chain has converged and the required number of samples has been collected.  Table SI.
The ABC and MD methods share several common features. Both of them rely on the simulation of epidemics using the current parameters, and on the calculation of a distance between simulated epidemics and the observed data. However, the ABC method considers any distance below the cutoff ǫ, while the MD method seeks to minimise the distance over a given number of simulations R. The iterative procedure used in MD method allows the use of an additional parameter analogous to ǫ in ABC method to be avoided. In both cases, the exact posterior is recovered in the proper limit (ǫ → 0 and R → ∞, respectively).

C. Markov Chain Monte-Carlo (MCMC) method with data augmentation
In order to implement the data-augmented MCMC method, it is necessary to calculate the explicit form of the likelihood È(D|θ) (with θ = (T, τ 0 , k)). We sketch here the main steps of the calculation.
Let I be the set of the dots that are colonised before the end of the experiment (at time t end = 21 days), and U be the set of those that are still uncolonised at t = t end . Assume first that we know the times of colonisation t j of each dot j ∈ I. The data then consist of the vector t of colonisation times, plus the set of uncolonised dots, i.e. D = (t, U ). The nearest neighbours of dot j form the set N j , and the potential donors of j form the subset S j ⊆ N j . If j ∈ U , then S j contains the colonised neighbouring dots, i.e.
S j = {i : i ∈ N j ∩ I}. If in contrast j ∈ I, S j contains the neighbouring dots that are colonised before Given these definitions, the likelihood function can be written as the product of the contributions from individual dots j: where f I j (t j ), is the p.d.f. for the colonisation times t j , and P U j (t end ) is the probability for dot j to be uncolonised by the end of the experiment. These contributions can be calculated explicitly as follows.
The probability that a dot j ∈ I has not been colonised by a given neighbour i ∈ S j by time t j is given by the survival function S(t j − t i ) (see Eq. (S.2)). Hence, the probability P I j (t j ) that dot j is still uncolonised at time t j is given by the product over all i ∈ S j : where we used the relation S(t) = exp(− t 0 φ(u)du). The p.d.f. for t j is then given by: Likewise, a dot j ∈ U is still uncolonised by time t end when transmission did not occur from any of its neighbours i ∈ S j , yielding the probability: where the integral is performed over the (high-dimensional and, in general, very complex) space T(d) compatible with the observed data.
Since the high-dimensional integrals introduced above are analytically intractable, numerical techniques are commonly used to sample values from the posterior distribution π(θ|D). The MCMC method consists in implementing a Markov chain for θ that has π(θ|D) as stationary distribution. A large literature exists on this subject (see e.g. Ref. [5]), to which the reader is referred for details. In our case, a Metropolis-Hastings algorithm has been used to build the Markov chain.
The algorithm can be summarised as follows: MCMC.1 Set n = 0 and choose the initial value θ 0 of the parameter vector.
MCMC.2 Generate a candidate value of the vector θ ′ from a proposal distribution q(θ ′ |θ n ).
MCMC.5 Set n = n + 1 and return to MCMC.2 until the chain has converged and the required number of samples has been collected.
In order to deal with the unobserved colonisation times t, we used data augmentation [2]. Their values were treated as parameters to estimate, i.e. the original parameter vector θ was expanded (augmented) to (θ, t).
New colonisation times were then proposed and accepted/rejected within the same Metropolis-Hastings algorithm. Additional care had to be taken in this case, since at each step a pathway of transmission must exist between the dot inoculated at time t = 0 and all the other colonised dots in the system (see the discussion in [2] and [6]).
Independent priors were used for all the parameters, so that π(θ) = π(T )π(τ 0 )π(k). A noninformative uniform prior was used for the transmissibility, i.e. π(T ) = U (0, 1). For the other two parameters, exponentially distributed priors were used, i.e. π(τ 0 ) = Exp(1), π(k) = Exp(1). Such choice was made in order to exploit our prior knowledge, i.e., respectively, that the typical time scale of the fungal spread in our system is of the order of days, and that the waiting-time distribution is not step-like (i.e., k is not too large; note that such assumption is opposite to that of the RF model). We checked that the posterior distribution was not too sensitive to changes up to a factor 2 in the parameters of the exponential priors. A uniform prior was used for each augmented colonisation time.
Every chain was run for 10 5 MCMC steps discarding an initial burn-in period of 10 3 steps. We

III. ADDITIONAL COMPARISONS OF FITTED MODELS WITH EXPERIMENTAL DATA
In the main text, we presented a test of the goodness of fit of models to data based on the mean squared distances d 2 c and d 2 f . The purpose of this section is to give a more visual comparison between fitted models and observations based on the cumulative incidence.
For every model used, we obtain the statistics for the incidence corresponding to fitted models by sampling the parameters θ from the joint (exact or approximate) posterior, π(θ|D). This procedure gives a p.d.f. ρ(C|t) for the incidence C at any given time t. The dispersion of ρ(C|t) is associated both with the stochasticity in the simulated model for each value of the parameters, θ, and the dispersion for the values of these parameters given by π(θ|D). Strictly speaking, the comparison of the experimental incidence with that obtained by methods A-D based on the RF model makes sense only when the stochastic nature of the process is taken into account. Indeed, the fine details of the two types of processes are different: while the experimental curve C(t) corresponds to a discrete sampling of a continuous-time process (discrete set of observations), ρ(C|t) gives the statistics of effective discretetime processes with random values of transmissibility corresponding to θ being drawn from π(θ|D). As can be seen from the figures, the p.d.f ρ(C|t) obtained by methods C, D, and E give a better description of the observed data for most of the replicates. On the other hand, the p.d.f. obtained by methodology F are often not able to capture the global trend of the incidence (see, e.g., replicate 3). This suggests that the focus of methodology F on finer details of the evolution may prevent the predictability of the invasive properties of the system. As discussed in the main text, this might be due to the negative interplay between the simplicity of the CT model and the individual-based description of the data-augmented MCMC method for inference.  Table 1 of the main text). In all the panels, the ridge (bold solid line) corresponds to the median of ρ(C|t). The grey-scale shaded areas are the 20% (darker), 40%, 60%, and 80% (lighter) percentiles around the median.

IV. FORECAST OF THE INCIDENCE
In the main text, we quantified the differences between predictions and observations for all methods and fungal invasion experiments in terms of the quantities ∆c and ∆F (cf. Fig. 6 of the main text).
In this section, we give a more visual comparison between observations and predictions based on the temporal incidence, C(t). Figs. S4 and S5 show the comparison for the six replicates available in the experiments with a = 10 mm and a = 12 mm, respectively. Predictions of the incidence between days 11 and 21 are based on estimates of the transmission parameters made during the first 10 days with methodology C. As can be seen, the estimated incidence provides a reasonable statistical description for the observed incidence for all the replicates. The quantity ∆c quantifies the rms distance between the observed incidence and the predicted incidence obeying the p.d.f ρ(C|t) (given by the grey-scale shaded area in Figs. S4 and S5).

A. Homogeneous transmission of infection
Here, we give numerical support to the claim made in the main text that the most probable estimate for transmissibility, T ⋆ , corresponding to the maximum of the probability density function (p.d.f.) ρ( T ), does not differ significantly from the actual transmissibility, T . This is shown in Fig. 6 where the estimate for the transmissibility, T , is plotted as a function of the transmissibility T . The estimates have been obtained for many different SIR numerical epidemics (∼ 10 4 ) with T ∈ [0, 1]. The analysis has been restricted to epidemics with final size N R (i.e. the number of removed hosts) greater than a certain cut-off, N 0 , in order to avoid estimates for small epidemics giving a poor estimator for transmissibility.
Excluding small epidemics from the analysis also makes sense from a practical point of view because they are not a threat in terms of invasion. We have checked that the statistics for T do not depend The first moment of ρ( T |T ) e gives an estimate for the mean T ⋆ e averaged over stochastic realisations.
The dependence of T ⋆ e on T is shown by the continuous line in Fig. S6. The dispersion of ρ( T |T ) e , shown by the shaded region in Fig. S6 as a function of T , contains contributions from both the width of each individual distribution, ρ i ( T |T ), and dispersion of the maxima for different replicates. In particular, the standard deviation, σ(T ), of ρ( T |T ) e is given by is the average over stochastic realisations of the variance σ 2 i of ρ i ( T |T ). The quantity σ 2 ⋆ is the variance of T ⋆,i over stochastic realisations, calculated as As can be seen from Fig. S6 the actual value for the transmissibility is statistically well described by the distribution ρ( T |T ) e . The mean for the most probable estimate for the transmissibility, T ⋆ e , is in good agreement with the actual transmissibility. In particular, T ⋆ e provides an excellent estimate for

B. Heterogeneous transmission of infection
In the main text, we have analysed the epidemics in model systems with homogeneous transmissibility for all pairs of connected hosts. Realistic populations of hosts exhibit inherent heterogeneity in transmissibility and it is crucial to understand its effect on the prediction method introduced in the main text.
In order to study the predictability of invasion for epidemics with heterogeneity in transmission of infection, we consider a simple but generic situation in which transmission is heterogeneous due to variability in the infectivity, I, and susceptibility, S, of hosts. As a first approximation, the rate of infection from an infected donor host with infectivity I d to a susceptible recipient host with susceptibility S r is defined as β d-r = I d S r [8]. We assume that I d and S r are independent random variables distributed according to truncated normal distributions, N (Ī, σ 2 I ) and N (S, σ 2 S ), respectively, which are the same for each host [9]. The mean values,Ī andS, provide an effective measure of the mean strength of the transmissibility while the standard deviations, σ I and σ S , characterize the degree of heterogeneity. The multiplicative form of the infection transmission rate β d-r = I d S r brings correlations in transmissibilities, [10]. Indeed, all the transmissibilities from a donor are affected by the value of I d and thus they are not independent. Similarly, all the transmissibilities to a recipient are influenced by its susceptibility S r and thus also correlated. Such correlations make the invasion probability for heterogeneous system to be dependent on the whole set of transmissibilities {T d-r }. In spite of that, the method based on a single effective transmissibility, T , and Eq. [1] in the main text is still applicable and useful. In realistic situations, the transmissibility is often assumed to be homogeneous because the precise degree of heterogeneity in transmission of infection is unknown and difficult to infer in detail.
In order to test our methodology in heterogeneous systems with different strengths of transmissibility,  each of the numerical epidemics in a system of size L = 51 (see Fig. 2(a) in the main text) by using Eq. [1] in the main text, the estimated ρ( T ), and the probability of invasion P hom ( T ; L) for systems with homogeneous transmissibility equal to T . Note that here P hom ( T ; L) corresponds to the function denoted as P inv ( T ; L) in the main text. The notations have been changed in order to distinguish between the probability of invasion in homogeneous systems and the probability of invasion in the presence of heterogeneity, denoted as P het ( T ; L). The results of our estimations are shown in Fig. S8.
The relation between P inv (L) and P hom ( T ; L) is very similar to the relation reported in the main text for epidemics with homogeneous transmission. Indeed, for epidemics with low transmissibility, P inv (L) typically overestimates P inv . In contrast, for more invasive epidemics, P inv (L) underestimates P inv for most of the possible effective transmissibilities. Although this comparison has some interest, in the current situation it makes more sense to compare the estimated probability of invasion with the actual probability of invasion in heterogeneous system, P het ( T ; L), that can be calculated numerically and shown by the line with squares in Fig. S8. The comparison of this line with the shaded region reveals that the estimated P inv (L) overestimates the actual probability of invasion in most of the situations both for cases with low and high transmissibility. Therefore, the estimates of P inv typically provide safe bounds for the probability of invasion. In fact, the larger the heterogeneity in susceptibility and/or infectivity, the safer the bound is. This is a consequence of the inequality P inv ( T ) ≥ P het ( T ) that holds for any given value of T under quite general conditions due to the existence of correlations in transmission induced by heterogeneity in the transmission rates [8,11]. Indeed, Fig. S8 shows that the inequality holds for the numerical experiments considered here with T = T (i.e. the line corresponding to the heterogeneous case marked by the squares is below the line marked by the circles for homogeneous system). These results are particularly encouraging for analysis of realistic epidemics in which a certain degree of heterogeneity in susceptibility and infectivity of hosts is expected to be ubiquitous.

VI. DIFFERENCES BETWEEN P inv ANDP inv
In this section, we discuss the origin of the difference between the estimated probability of invasion, P inv (L), evaluated at the most probable transmissibility, T ⋆ , and the actual probability of invasion P inv ( T ; L) one would obtain if the transmissibility of the epidemic was known exactly.
To start with, recall that the relation betweenP inv (L) and P inv ( T ; L) is given by Eq. [1] in the main text. As we have seen, ρ( T ) is a peak-shaped function approximately symmetric about the peak position at T ⋆ . Therefore, the peak region will mainly contribute to the integral in Eq. is sufficient to estimate the deviation of P inv from P inv . Indeed, substitution of Eq. (S.10) into Eq. [1] of the main text gives: The distribution of T is approximately symmetric around the maximum, i.e. ρ( T − T ⋆ ) ≃ ρ( T ⋆ − T ) (see Fig. 2(b) in the main text), and thus the term containing P ′ inv in Eq. (S.11) is negligible in comparison with other terms in the sum. This means that P inv (L) > P inv ( T ⋆ ; L) if P ′′ inv ( T ⋆ ; L) > 0 P inv (L) ≤ P inv ( T ⋆ ; L) if P ′′ inv ( T ⋆ ; L) ≤ 0 , (S.12) where we have taken into account that the last integral in Eq. (S.11) is positive. The above inequalities demonstrate that the value and sign of P inv (L) − P inv ( T ⋆ ; L) depends on the curvature of P inv ( T ; L) around T = T ⋆ which is given by P ′′ inv .