## Abstract

Using digitized images of the three-dimensional, branching structures for root systems of bean seedlings, together with analytical and numerical methods that map a common susceptible–infected–recovered (‘SIR’) epidemiological model onto the bond percolation problem, we show how the spatially correlated branching structures of plant roots affect transmission efficiencies, and hence the invasion criterion, for a soil-borne pathogen as it spreads through ensembles of morphologically complex hosts. We conclude that the inherent heterogeneities in transmissibilities arising from correlations in the degrees of overlap between neighbouring plants render a population of root systems less susceptible to epidemic invasion than a corresponding homogeneous system. Several components of morphological complexity are analysed that contribute to disorder and heterogeneities in the transmissibility of infection. Anisotropy in root shape is shown to increase resilience to epidemic invasion, while increasing the degree of branching enhances the spread of epidemics in the population of roots. Some extension of the methods for other epidemiological systems are discussed.

## 1. Introduction

The identification of efficient strategies is of great practical importance in controlling the spread of epidemics in populations of plants and animals and in social and computer networks (Nâsell 2002; Pastor-Satorras & Vespignani 2003, 2004). Many previous studies from different disciplines, such as epidemiology, mathematics, physics, chemistry and social sciences, have attempted to address this problem using experimental, numerical and analytical approaches (e.g. Liggett 1985; Marro & Dickman 1999; Diekmann & Heesterbeek 2000; Murray 2002; Keeling & Rohani 2007; Jackson 2008).

Numerous mathematical models exist for the description of epidemics in populations (e.g. Nâsell 2002 and references therein). Typically, a population is represented by a network in which nodes are associated with hosts and links between nodes describe possible transmission paths for diseases. The hosts can be in several states such as susceptible (S), infected (I) and recovered (R) (or equivalently removed). In two prototype epidemiological models, SIS and SIR, the hosts change their state according to certain rules. For example, in the SIS model, a susceptible host (S) can be infected by its neighbour (I) and then again become susceptible, thus following the sequence S → I → S. The SIR model differs from the SIS only in the last step of the sequence, S → I → R, at which an infected host is effectively removed from the epidemic either by death or by becoming immune to further infection. Extensions and modifications of these two models have been proposed to describe various features of realistic epidemics (e.g. Dybiec *et al.* 2004).

Both models follow stochastic dynamical rules meaning that either one or both of the dynamical processes, infection and/or recovery, occur randomly, thus reflecting the inherent complex nature of the hosts and communications between them. Depending on the relative value of the infection and recovery rates, the epidemic can be in different regimes: non-active (epidemic does not cover the major part of the population) and active (outbreak of epidemic). The transition between these two regimes occurs at a critical point characterized by the critical value of a control parameter (e.g. the ratio of recovery and infection rates) and can be described as a continuous phase transition belonging to the universality class of directed or isotropic percolation for SIS and SIR models, respectively (Marro & Dickman 1999; Hinrichsen 2000; Ódor 2004). Universality at the threshold refers to universal critical exponents describing scaling behaviour of, for example, correlation length in spatial fluctuations of hosts in a certain state. In contrast, the threshold value of the control parameter is system dependent and can vary significantly between models belonging to the same universality class.

The location of the threshold in the parameter space and the evaluation of the probability of invasion above the threshold provide information about the vulnerability of the system to epidemic outbreaks. The value for the threshold is known analytically only for some simple models: for example, the threshold is known for both the SIR and SIS models in the case of networks with homogeneous transmission rates and given degree distributions (Pastor-Satorras & Vespignani 2003). More usually, however, numerical methods are required to identify both the probability and threshold for invasion. The derivation of analytical insight is further challenged for many realistic systems that are characterized by inherently heterogeneous parameters, which may also be correlated in space.

Heterogeneity can occur among the nodes and connections in epidemiological networks. Thus individual host characteristics, such as host type and age, recovery time, infectivity and immunity may differ; there may also be variation in the transmissibility of infection between pairs of infected and susceptible hosts. Recently, Pérez-Reche *et al.* (2010) considered heterogeneity among nodes by analysing the influence of host morphology on the properties of epidemics in a set of two-dimensional retinal ganglion cells placed on the nodes of a regular two-dimensional lattice. They showed that heterogeneity associated with host morphology significantly influenced the invasion threshold, making a system of morphologically complex hosts less vulnerable to epidemic invasion (Pérez-Reche *et al.* 2010). Pérez-Reche *et al.* (2010) also showed that, although the set of neural cells exhibited heterogeneity, there was little evidence for correlation in the transmission of infection between pairs of contiguous hosts. Transmission could therefore be successfully described (Pérez-Reche *et al.* 2010) by an equivalent effective mean-field homogeneous system (Sander *et al.* 2002). Many systems, however, exhibit greater heterogeneity than the two-dimensional ganglion system, especially in the degree of anisotropy and in the transmissibility between hosts. In this paper, we advance the analyses to three-dimensional morphologically complex and anisotropic hosts (plant root systems) and demonstrate that high anisotropy in three-dimensional root systems can bring non-negligible correlations in pair-wise transmission of pathogens from infected to susceptible root systems. Such correlations lead to a breakdown in the mean-field description and make these systems safer for epidemic outbreaks.

Below, we focus on the spread of soil-borne disease through populations of plants in which the transmission of infection occurs between contiguous plants and is mediated by the complex morphology of branching root structures (Gilligan *et al.* 1994). Specifically, the aims of the paper are (i) to study analytically and numerically the spread of epidemics in a population of plants with realistically complex three-dimensional morphology, (ii) to locate the associated threshold for invasion and subsequent epidemic spread, (iii) to predict how the invasion threshold depends upon the shape of the three-dimensional hosts, and (iv) thus to suggest how to control the vulnerability of the system to pathogen invasion by changing the host morphology. In this paper, we consider a population of plants, each comprising an ensemble of roots, with the plants arranged on a regular two-dimensional lattice to mimic the spatial arrangement of plants in a field. For concreteness, we deal with an SIR model in which infected plants eventually die (i.e. permanently go to the R-class). The SIR model can be applied to a wide range of soil-borne pathogens (Hornby 1990) spreading through populations of roots (Gilligan 1990, 2002). In what follows, we analyse how the morphology of the root systems of bean (*Phaseolus vulgaris*) seedlings affects the probability of epidemic invasion of a common and ubiquitous soil-borne, fungal pathogen, *Rhizoctonia solani*, which spreads by mycelial extension from infected to susceptible hosts (Bailey *et al.* 2000).

Our main findings are the following. (i) The morphology of the hosts represented by three-dimensional roots is an important factor for the SIR epidemic. Two morphological features, namely anisotropy and disorder, are of particular importance for the spread of epidemics. Anisotropy refers mainly to the deviation of the major root stems from the vertical direction. Disorder is a characteristic related to the degree of root branching and in particular to the number of secondary branches in the root. (ii) Anisotropy brings short-range correlations to the disease transmission that makes the system less vulnerable to epidemic invasion. (iii) Disorder diminishes the effect of anisotropy and thus increases the vulnerability of the system to an epidemic outbreak.

The structure of the paper is as follows. In §2, the SIR model is introduced for the root system placed on a square lattice. The main results for invasion probability and analysis of correlations in transmissibilities are presented in §3. Discussion and conclusions are given in §4.

## 2. Model

### 2.1. SIR process

In the SIR model used below, the hosts can be in one of the three states, susceptible (S-class), infected (I-class) or recovered and *permanently* immune to infection (R-class). Here, a host refers to a plant root system. Each plant is located on the node of a square (for concreteness) lattice with links between nearest neighbours only. An infected host remains in such a state for a fixed time *τ* and then ceases to be infectious (i.e. permanently goes to the R-class). The recovery time is chosen to be the timescale for the problem and set to *τ* = 1. During the infectious period, while times are in the range *t* ≲ *τ*, an infected node can stochastically, by a Poisson process, transmit the pathogen and hence cause disease to its nearest neighbour with a rate *β*. Under these rules, the probability, *T*, that disease is passed from an infected host to its susceptible neighbour before that host moves to the R state is given by the following expression (Grassberger 1983):
2.1

This value, called the transmissibility, quantifies the transmission process between neighbouring hosts. At a large scale, *T* determines the probability *P*_{inv} that a large outbreak occurs, with the pathogen eventually invading the system (invasion probability). The probability *P*_{inv} is difficult to calculate analytically, and typically it is computed numerically by generating many stochastic realizations of epidemics and counting the fraction of those that span the system across all borders.

In the language of phase transitions, *T* is the control parameter for the homogeneous SIR model (it is varied to study the transition) and *P*_{inv} is the order parameter. If transmissibility is less than some critical value, *T* < *T*_{c}, then epidemic invasion cannot occur in an infinite system and the invasion probability for the epidemic is zero, *P*_{inv} = 0. If *T* > *T*_{c}, the invasion probability is finite, *P*_{inv} > 0, approaching unity when *T* → 1. The critical value of transmissibility coincides with the critical bond probability, *p*_{c}, for the bond percolation problem because the final cluster of R sites in the SIR problem coincides with one of the clusters of sites connected by occupied bonds in the bond percolation problem, if the bond probability is equal to the transmissibility (Grassberger 1983). In particular, *T*_{c} = *p*_{c} = 1/2 for a square lattice with homogeneous transmissibilities (Isichenko 1992; Stauffer & Aharony 1992).

An epidemic with heterogeneous uncorrelated transmission rates (for fixed value of *τ*) can also be mapped onto the bond percolation problem, so that, at criticality, the following condition should hold (Newman 2002; Sander *et al.* 2002; Kenah & Robins 2007; Miller 2007):
2.2
where *ρ*(*β*) is the probability density function for the transmission rate *β* and *p*_{c} is the bond percolation threshold for the lattice on which the SIR model is defined. Moreover, this mean-field relation is valid for the whole range of transmissibilities and the invasion probability in a heterogeneous system can be shown to be (Sander *et al.* 2002)
2.3
i.e. an epidemic in an uncorrelated heterogeneous system with transmissibilities {*T*_{ij}} is equivalent to an epidemic in a homogeneous system in which all the transmissibilities are replaced by the mean value of transmissibility 〈*T*〉.

The above statement is correct only in the case of independent transmissibilities between different pairs of hosts. Below, we use Monte Carlo simulations of multiple possible invasion events to demonstrate that for realistic hosts with complex morphology this assumption does not necessarily hold and a simple mean-field description of the SIR process can fail. Namely, the inherent short-range correlations in transmissibilities between different pairs of hosts are shown to change the threshold value for transmissibility.

### 2.2. Root systems as particular hosts

Realistic bean root systems were used as hosts in our model. Nine bean seeds, *N*_{r} = 9, were allowed to germinate for 2 days on cotton, and then transferred to an aquarium to grow hydroponically over 4 days, before being set into paraffin blocks which were then sliced into 0.1 mm thick layers. The layers were digitally scanned with pixels of a linear size equal to 0.1 mm, to produce three-dimensional images of the primary (tap) root and first-order lateral roots; hereafter referred to as *root systems* (see a typical image in figure 1*b*).

The root systems thus obtained (figure 1) were used as hosts in a lattice model to study the invasion of a soil-borne SIR epidemic spreading through a population of plants by mycelial growth from root system to root system. The lattice model was created in the following manner. An arbitrary root system (e.g. A in figure 1*a*) was selected from the nine experimentally available systems (see first three rows in figure 1*a*). This root system was rotated by a random angle 0 < *φ* < 2*π* about a vertical axis passing through the original bean and then placed on a node (with the bean-seed position coinciding with the node) *i* (*i* = 1, 2, … , *L*^{2}) of a square lattice of size *L* × *L* with lattice spacing *a* (e.g. root A′ in figure 1*a*). This process was repeated for all the nodes of the lattice. Several realizations (approx. 10^{3}) of so created lattice systems form a statistical ensemble, which has been analysed below.

A soil-borne epidemic can be modelled by introducing a pathogen, in the form of a small fungal colony typical of *R. solani* onto a root system near the centre of the lattice. The pathogen is assumed to infect the initial plant and to grow along and around the root system by mycelial extension exploring the soil nearby. If another root system is situated nearby, the fungal colony can reach it and eventually infect it. In this way, the fungal colony propagates through a population of roots.

The transmission rate between plants in the SIR process, which is effected by mycelial spread between root systems, is proportional to the effective overlap between different three-dimensional root systems. In order to characterize this overlap quantitatively, we first define the root number density (number of voxels per volume) in the following way. The roots (comprising the primary root and first-order lateral roots) in our model are represented by three-dimensional digital images. Each image, *i*, consists of a set of *N*_{i} voxels (three-dimensional pixels) with volume *V*_{voxel} = 10^{−3} mm^{3}, the locations of which are given by position vectors *r*_{ni} = (*x*_{ni}, *y*_{ni}, *z*_{ni}) with *n* = 1, … , *N*_{i}. The root number density that the image represents, *ρ*_{scan,i}(** r**), is equal to 1/

*V*

_{voxel}if

**coincides with the position of one of the voxels and zero otherwise. On the scale of the system, the number density of voxels can be written as a sum of Dirac delta functions, 2.4**

*r*Bearing in mind that the density should describe some additional volume around each root across which the fungus can spread though soil, known as the pathozone (Gilligan 1995), and which therefore mediates the transmission of infection from a donor (i.e. infected) to a recipient (susceptible) root system, we apply broadening to the images replacing the Dirac delta functions in equation (2.4) (or voxels) by their Gaussian representation, i.e.
2.5
with the broadening parameter *σ* = 1.5 (the value of *σ* is much smaller than the lattice spacing, i.e. *σ* ≪ *a*, so that the overlap with next-nearest neighbours and thus transmission of the pathogen to them can be ignored in the model) of the linear size of the voxel (px). This mimics the diffusive spread of a pathogen into the medium surrounding the root (Gilligan 1995).

The fungal colony can efficiently pass between different roots through the volume of soil if branches of both roots are present in that volume, i.e. if the roots overlap. Therefore, we assume that transmission between nearest-neighbour hosts *i* and *j* separated by unit cell vector takes place at a rate *β*_{ij}, proportional to the three-dimensional overlap *J*_{ij} between the two hosts, i.e
2.6
where
2.7

The value of *J*_{ij} (dimensionless) is the number of voxels of root *i* in the region of the overlap with root *j*. From the definition, *J*_{ij} = *J*_{ji}. The coefficient *k*, the infection efficiency, has the meaning of transmission rate per overlapping voxel and is set to be the same for all overlapping roots. This efficiency represents the ability of the pathogen to use a point of contact between two hosts to spread between them. It therefore takes into account all details of the transmission process not associated with the host morphology and specific to the particular host–pathogen system, such as susceptibility and infectivity, together with other factors such as environmental conditions that, in turn, influence transmission of infection. We therefore investigate the role of morphology in deciding the overall rate between two root systems, with all other properties described by variable values of *k*. Therefore, both the transmission rate, *β*_{ij}, between an arbitrary pair of the nearest-neighbour hosts and thus transmissibility,
2.8
(see equation (2.1)), depend on the transmission efficiency, *k*, being a control parameter, and the random overlap, *J*_{ij}, which depends on distance between hosts *a*, being another control parameter. The randomness in the overlaps is caused by the complex morphology of hosts. With these assumptions, all the information about a given spatial configuration of the system (a particular realization of roots and angles at each node) is contained in the set of overlaps {*J*_{ij}}, while all the information about the epidemiological properties of the system is contained in the set of transmissibilities {*T*_{ij}}.

Typical probability densities *ρ*(*J*_{ij}) and *ρ*(*T*_{ij}) for the ensemble of roots are shown in figures 2 and 3. It follows from figure 2 that *ρ*(*J*_{ij}) is a monotonically decaying function. The value of *J* roughly corresponds to the number of pixels in the overlap region. In the region of very small overlaps *J*_{ij} ≲ 1 (i.e. the overlap is either less than or of the order of unity) caused by Gaussian broadening effects for pixels,
2.9
with *α* ≃ 0.9 (see the inset in figure 2). For the most interesting range of 1 ≲ *J*_{ij} ≲ 10^{3}, the probability density decays according to equation (2.9) with *α* ≃ 0.5 (figure 2). For larger values of *J*_{ij}, the function *ρ*(*J*_{ij}) exhibits even faster decay.

The shape of the probability density for transmissibilities follows from the mapping given by equation (2.8). The semi-infinite range of *J*_{ij} is mapped onto a finite interval for transmissibilities, 0 ≤ *T*_{ij} ≤ 1, in such a way that all the large values of *J*_{ij} are accumulated around *T*_{ij} ≃ 1 and all small values of *J*_{ij} are concentrated around *T*_{ij} ≃ 0. As a consequence, the monotonic distribution *ρ*(*J*_{ij}) gives rise to a distribution with two maxima (figure 3). The relative weights of the peaks depend upon the value of the transmission efficiency, *k*. For small values of *k* ≲ *k*_{c}, the peak around zero is dominant (see the solid line in figure 3), ensuring the non-active regime for epidemics, i.e. there is no invasion. In contrast, for large values of *k* ∼ *k*_{c}, the peak around unity becomes dominant for the distribution of transmissibilities and the SIR process is active (see the dot-dashed line in figure 3). At criticality, *k* = *k*_{c}, the distribution of transmissibilities is approximately symmetric around *T* = 1/2 (see the dashed line in figure 3).

## 3. Results

In this section, we present our main results on the behaviour of the SIR epidemic for the ensemble of realistic three-dimensional root systems of bean plants placed on a regular lattice. We start with an evaluation of the invasion probability (order parameter) as a function of two control parameters, transmission efficiency and lattice spacing, paying special attention to the threshold values of the control parameters separating invasive and non-invasive regimes of epidemics. After that we compare the behaviour of the epidemic in our model with a mean-field system and find a significant role of correlations in transmissibilities of pairs of bonds separated by a relatively short distance. Finally, several ‘toy’ models are investigated in order to clarify the effect of correlations.

### 3.1. Phase diagram

One of the most important characteristics for epidemic spread is the probability of invasion. In the language of continuous phase transitions, the invasion probability, *P*_{inv}, is the order parameter defining the status of the epidemic to be in either the non-invasive (*P*_{inv} = 0) or invasive (*P*_{inv} > 0) regimes. In our model, the invasion probability depends on two control parameters—the transmission efficiency (specific to the host–pathogen system and independent of the spatial configuration) and the lattice spacing (affecting the set of overlaps {*J*_{ij}}), *P*_{inv} = *P*_{inv}(*k*, *a*).

In order to calculate *P*_{inv}(*k*, *a*) for a particular set of transmissibilities {*T*_{ij}}, we use the mapping to bond percolation and create many stochastic realizations in which each bond *i* − *j* is occupied with probability *T*_{ij}. The probability of invasion is defined in an infinite system as the probability that any given site belongs to an infinite, or spanning, cluster of sites. A single such spanning cluster will always exist above the invasion threshold in an infinite system, and will never exist below it. Therefore, the invasion probability is equal to the fraction of sites out of the whole system that belong to this cluster, when it exists, and zero otherwise. Since it is clear that there will not be an infinite cluster in a finite system (which is all we can model), we define a spanning cluster here as one which spans the lattice touching all four boundaries. Then for each stochastic realization of the bonds, we determine whether the spanning cluster exists, and what fraction of the system it makes up, thus obtaining an estimate of *P*_{inv} from that realization. The average of these estimates over all stochastic realizations gives the probability of invasion for the particular values of {*T*_{ij}}.

A different realization of the set {*T*_{ij}} would produce a slightly different value of the invasion probability. Therefore, we have averaged the invasion probability over different realizations of {*T*_{ij}}, to obtain the *configurational average* of the invasion probability 〈*P*_{inv}(*k*, *a*)〉. This value is representative because the invasion probability is a self-averaging quantity (Binder & Young 1986), as we have checked numerically.

The results of numerical simulations for the average invasion probability, 〈*P*_{inv}(*k*, *a*;*L*)〉, for a finite size, *L* × *L*, lattice are presented in figure 4, with *L* = 200. As expected, for large values of lattice spacing and small transmission efficiency, the system is unsuited to pathogen invasion. In contrast, small lattice spacings and high values of *k* enhance the probability of invasion and hence of an epidemic.

The invasion probability surface 〈*P*_{inv}(*k*, *a*;*L*)〉 shown in figure 4 depends on the lattice size *L*. When the system size becomes very large, *L* → ∞, the values of the probability of invasion tend to a certain limiting surface. The line of the critical points *a*_{c} = *a*_{c}(*k*_{c}) (the solid line in the horizontal plane in figure 4) separating non-invasive (lim_{L→∞}〈*P*_{inv}(*k*, *a*|*L*)〉 = 0) and invasive (lim_{L→∞} 〈*P*_{inv}(*k*, *a*|*L*)〉 > 0) regimes is of particular importance because it allows the identification of the safe (non-invasive) parameter region.

The data points on the phase boundary *a*_{c} = *a*_{c}(*k*_{c}) were obtained by finite-size scaling collapse in the following way. The dependence of the order parameter on the system size around the critical point (e.g. *k*_{c}) for continuous phase transition is well established (Isichenko 1992),
3.1
where *P̃*_{inv} is the scaling function and *β* and *ν* are the universal scaling exponents. The scaling function can be found by scaling collapse of 〈*P*_{inv}(*k*, *a*;*L*)〉 for several system sizes using *k*_{c} and scaling exponents as fitting parameters (Heermann & Stauffer 1980). Technically this is achieved by plotting 〈*P*_{inv}(*k*, *a*;*L*)〉*L*^{β/ν} versus (*k* − *k*_{c})*L*^{1/ν} for several values of *L*, and finding the values of *ν*, *β* and *k*_{c} that lead to a collapse of all the curves for different *L* to a single one, which corresponds to *P̃*_{inv}. Such a finite-size scaling analysis was applied to evaluate the critical values of *k*_{c} for given lattice spacings and the results are shown by the line in the horizontal plane in figure 4. The scaling collapse also gives the values of the scaling exponents, *β* = 0.13 ± 0.008 and *ν* = 1.32 ± 0.05, which coincide with those for the isotropic percolation universality class (Isichenko 1992). We have also performed the scaling collapse for the susceptibility (Heermann & Stauffer 1980) and found the expected values for the corresponding scaling exponents and the invasion threshold. Therefore, we have confirmed that the continuous phase transition for invasion probability in the system of roots (with short-range correlations discussed below), as expected, is similar to that in isotropic percolation, i.e. it belongs to the isotropic percolation universality class.

### 3.2. Correlations

The invasion probability presented in figure 4 for a heterogeneous lattice with transmissibilities {*T*_{ij}} has been calculated numerically by estimating the mean fraction of sites belonging to the spanning cluster in the bond percolation problem. Alternatively, we could use a mean-field assumption (equation (2.3)), which suggests that the heterogeneous system is equivalent to a homogeneous one with average transmissibility 〈*T*〉, if there are no correlations between transmissibilities. This would prove that the spread of the soil-borne pathogen in the root system is not being affected by the heterogeneity of the hosts. To test the validity of the mean-field approach for our model, in figure 5, we show the dependence of the invasion probability on mean transmissibility both for heterogeneous (solid curve) and for homogeneous mean-field (dashed curve) systems. Evidently, the critical value of the transmissibility in the real-root system is greater than that in the mean-field system, and thus the real-root system with inherent heterogeneity in transmissibilities is less vulnerable to epidemic invasion than the homogeneous system with the same MEAN transmissibility. This implies the presence of correlations in transmissibilities between hosts in the system (Sander *et al.* 2002; Pérez-Reche *et al.* 2010). Note that these correlations are of a different nature from spatial–temporal correlations in the density of infected/removed hosts observed during the propagation of disease. The short-range correlations in transmissibilities make an epidemic less likely in the heterogeneous system under consideration (solid curve in figure 5) than in a homogeneous mean-field system (dashed curve) for values of *P*_{inv} ≲ 0.95.

It should be noted that, for relatively large values of the mean transmissibility (e.g. 〈*T*〉 0.62 in figure 5), the invasion probability for a heterogeneous system is higher than that for a homogeneous one. This is a consequence of negative correlations in transmission (cf. figure 5). This behaviour is in contrast with that for systems with non-negative correlations in transmission, such as a system with heterogeneous recovery times. Indeed, in this case it can be rigorously shown that the probability of invasion in a homogeneous system cannot be smaller than that for heterogeneous systems for any value of 〈*T*〉 (Kuulasmaa 1982; Cox & Durrett 1988; Kenah & Robins 2007; Miller 2007). However, this theorem does not apply to systems with negative correlations, such as the root system studied here.

The difference in the invasion curves for heterogeneous and homogeneous systems (cf. solid and dashed curves in figure 5) is a clear indication of the correlations in transmissibilities for heterogeneous systems. Indeed, the spatial correlation function for transmissibilities,
3.2
plotted in the inset to figure 5, shows the presence of negative short-range (between nearest bonds) correlations (see the dip at *r*/*a* = 1 in the inset of figure 5). In equation (3.2), the distance *r* is taken between midpoints of different bonds *i* − *j* and *k* − *m* and the averaging is performed over all bond pairs separated by the same distance.

Therefore, the SIR process in the system of real roots can be mapped onto the bond percolation problem with short-range correlations between bond probabilities. In the next subsection, we reveal the origin of such correlations.

### 3.3. Disc-based models

In order to understand the origin of the short-range correlations in transmissibilities between roots, we study three simplified ‘toy’ models of differing complexity. In all these models, each root system *i* is replaced by a two-dimensional horizontal ‘disc’, i.e. by the number density *ρ*_{i}(** r**) with the centre of mass at

*r*_{i}. The centres of the discs are randomly displaced from the lattice nodes characterized by position vectors

*r*_{i}

^{(0)}, i.e.

*r*_{i}=

*r*_{i}

^{(0)}+

*Δ*

*r*_{i}, where the displacement vectors are normally distributed,

*Δ*

*r*_{i}∼

*N*(0,

*σ*

_{Δ}

^{2}), with zero expectation value and variance,

*σ*

_{Δ}

^{2}. This normal distribution of

*Δ*

*r*_{i}mimics the anisotropy in the shape of real-root systems, i.e. the fact that the centres of mass of roots are displaced with respect to the seed (node) positions. The toy models then differ among each other by the form of the function used for

*ρ*

_{i}(

**) and therefore by the value of overlaps between two discs.**

*r*In the first model (figure 6*a*), all the discs have the same radius, *R*, and constant number density, *ρ*_{i} (** r**) =

*ρ*

_{0}if |

**−**

*r*

*r*_{i}| ≤

*R*, and zero density otherwise. The overlap between two discs

*i*and

*j*with centres separated by distance

*d*

_{ij}= |

*r*_{j}−

*r*_{i}| reads 3.3 if

*d*

_{ij}/2

*R*≤ 1 and zero otherwise (with the px area

*V*

_{px}= 1). This is one of the simplest models which is equivalent to a representation of the real-root systems by identical uniform vertical cylinders projected onto the horizontal plane.

In the second model (figure 6*b*), the density is normally distributed around *r*_{i} with variance *σ*_{R}^{2}, i.e.
where *ρ̃* is a dimensionless normalization constant, and the overlap between two dispersed discs represented by Gaussian functions is
3.4
The second model accounts for the property that the number density in real roots is a decaying function with distance from the main tap root.

The third model (figure 6*c*) mimics the angular irregularity in the shape of real-root systems, i.e. the fact that overlap significantly depends on the angle *φ* at which the roots are placed at the lattice nodes (figure 1). In this model, we assume that the overlap has two contributions,
3.5
where *J*_{ij}^{(1)} is given by equation (3.3) as in the first model and a random normally distributed constituent, *δ**J*_{ij} ∼ *N*(0, *σ*_{J}^{2}), with the variance, *σ*_{J}^{2}, which can depend on *J*_{ij}^{(1)}.

In order to illustrate that irregularity in shape alone is insufficient to cause the observed effect, we consider a fourth model, which is a special case of the third model, where there is no displacement of the discs (i.e. we set *σ*_{Δ}^{2} = 0). Here we see no correlations (see the dot-double dashed line with left triangles in the inset in figure 7), and thus the system can be approximated by an equivalent homogeneous system (the dot-double dashed line with left triangles coincides with the solid line marked by crosses in figure 7).

The disc-based models depend on several parameters, such as *R*, *σ*_{Δ}^{2}, *σ*_{R}^{2} and *σ*_{J}^{2}. The two normalization constants *ρ*_{0} and *ρ̃* can be incorporated into the value of the infection efficiency *k* and are set, for convenience, to *ρ*_{0} = *ρ̃* = 1. The values of the disc radius, *R*, in the first model and the width of the Gaussians in the second model, *σ*_{R}, can be estimated from the mean value (among all *N*_{r} roots) of the variance of the real-root density projected onto the horizontal plane, i.e.
3.6
where
3.7
with *ρ*_{i} (** r**) given by equation (2.5) and suffix ⊥ referring to the vector component perpendicular to the vertical axis. This value is found to be

*R*=

*σ*

_{R}≃ 73 px

^{1/2}. Similarly, the value of

*σ*

_{Δ}

^{2}was estimated as 3.8

The value of *σ*_{J}^{2} can be estimated from the variance of the overlap *J*_{ij}(*φ*_{i}, *φ*_{j}) between two real roots, which stochastically depends on their orientations characterized by *φ*_{i} and *φ*_{j} (figure 1), and it is found to be *σ*_{J}/〈*J*_{ij}^{(1)}〉 ≃ 1.35.

Having the estimates of all the parameters for the toy models, it is straightforward to evaluate numerically the dependence of the invasion probability on mean transmissibility. The results presented in figure 7 clearly demonstrate that, for all the models exhibiting anisotropy, the critical point *T*_{c} is shifted to the right of that for the mean-field system and the invasion probability is also reduced in a rather wide range, 0 ≤ *P*_{inv} ≲ 0.95, when compared with the homogeneous system (note that models 2 and 3 reproduce the real-root system reasonably well). This is a consequence of the short-range correlations in transmissibilities present in the disc-based models. The bigger shift for the first model when compared with less significant shifts for models 2 and 3 reflects the respective decrease in the scale of correlations for the latter models (cf. dashed, dot-dashed and double dot-dashed curves in figure 7 and corresponding curves in the inset). The relatively large correlations in the first model are due to sharp disc boundaries so that the overlaps are significantly larger (smaller) for spatially close (distant) discs when compared with the mean value. This difference is diminished by the Gaussian broadening of the disc density so that the correlations for the second model are less pronounced. The angular fluctuations of density (model 3) result in a similar reduction of correlations found for the solid disc model.

We conclude from the analysis of the toy models that anisotropy in the root shape accounts for most of the short-range correlations in transmissibilities between the root systems. Anisotropy arises mainly from the horizontal shift of the centre of mass with respect to the lattice node (seed) location. To the greatest extent, this effect is picked up by the first model, which displays the most significant effect of correlations. Dispersion and irregularity in the root shapes mimicked by models 2 and 3, respectively, reduce the level of correlations. Finally, if irregular discs are not displaced from the lattice nodes (model 4), the correlations completely disappear and the model is equivalent to the mean-field one.

## 4. Discussion

The occurrence of well-defined phase transitions in determining whether or not a pathogen invades through a lattice of homogeneous hosts is well established (Grassberger 1983). Following initial work by Bailey *et al.* (2000) on the spread of the soil-borne fungus *R. solani*, critical percolation phenomena and phase transitions in lattice-based systems have been supported by experimental analysis. More recently, the concept has been extended to the transmission of soil-borne infections through lattices with missing sites (Otten *et al.* 2004) and through animal populations, exemplified by the spread of plague through gerbil populations in lattices of interconnecting burrows (Davis *et al.* 2008). Here, we have extended the theoretical analysis of critical percolation phenomena to take account of the spatial correlations that occur in the set of morphologically complex three-dimensional hosts represented by plant root systems.

We have used digitized images of the root systems of bean seedlings to provide realistic examples of a complex morphological branching structure to test the methodology. However, our results and methodological approaches, in particular the new disc-based models, provide a basis for further analysis of the transmission of infection and disease through a broad range of morphological structures.

We have shown that the invasion probability (*P*_{inv}) of a root-infecting pathogen spreading through a population of plants with morphologically complex, overlapping three-dimensional root systems depends upon the transmissibility of infection between neighbouring plants (equation (2.1)). We have resolved the transmission into two components, the transmission efficiency (*k*) between nearest neighbours and the lattice spacing between plants (*a*), which control the degree of overlap. The transmission efficiency is a measure of the probability of infection arising from overlap between an infected and a susceptible root system. Here, we assume a constant transmission rate per overlapping voxel for overlapping roots but this assumption could be relaxed. The extent of the overlap is broadened to allow for the ability of the pathogen to explore the soil immediately surrounding a root, known as the pathozone (Gilligan 1983). The spacing between plants controls the degree of overlap. We have used the mapping to bond percolation in order to calculate the invasion probability for a set of transmissibilities, showing evidence for a marked phase transition that depends upon *k* and *a* (figure 4). We conclude that the heterogeneities in transmissibilities arising from correlations in the degree of overlap between neighbouring plants render a population of root systems less vulnerable to epidemic invasion than a corresponding uncorrelated homogeneous system.

For concreteness, we have analysed the square-lattice topology only. Any other regular lattice topologies can be treated in a similar way. The invasion curve, *P*_{inv} (〈*T*〉), is specific for the lattice topology but the effect of correlations in transmissibilities between root systems will be qualitatively the same as that found for the square lattice. In modelling the spread of an SIR epidemic through heterogeneous, overlapping root systems, we have used several simplifications. One of these is related to the assumption of identical recovery times for all the hosts. If this assumption is relaxed and the recovery times are distributed according to a certain probability density function, then the invasion probability in such a system is lower than that in the system with homogeneous recovery times (Kuulasmaa 1982; Cox & Durrett 1988). We have demonstrated this using several realistic distributions for recovery times in root systems. The results presented in figure 8 clearly show a shift of invasion curves to the right, thus indicating that systems with heterogeneity in recovery times are indeed more resistant to SIR epidemics. The main reason for such increased resistance is that the heterogeneity in recovery times brings positive correlations to the system (see inset in figure 8; Kuulasmaa 1982; Cox & Durrett 1988; Kenah & Robins 2007; Miller 2007).

Another assumption in our approach is that the root anisotropy is in a random direction for each root. However, in a real planting, large-scale resource gradients can lead to a predominant orientation of all root systems. This results in anisotropy of the overlaps *J*_{ij}, meaning that the root overlaps along a certain direction are greater than those in other directions. In fact, the system becomes more ordered, typically exhibiting larger transmissibilities in the direction of resource gradients. For illustration, we considered an extreme case for which an epidemic spreads through a system of anisotropically identical roots placed in the same orientation (without rotation) on all the nodes of a square lattice. The results are presented as the solid curve marked by triangles in figure 8 from which it follows that the system of identical anisotropic roots is more vulnerable than the system of randomly rotated roots (cf. the dashed curve with the solid one), meaning that the external gradients may reduce the resilience of the system. This is an expected effect because the SIR epidemic in the system of identical anisotropic roots can be mapped onto anisotropic percolation characterized by two transmissibilities *T*_{x} ≠ *T*_{y} along different lattice directions. It is known (Sykes & Essam 1964) that the invasion threshold in such a system coincides with that in a homogeneous one, i.e. given by the equation *T*_{c} = 〈*T*〉 = (*T*_{x} + *T*_{y})/2 = 1/2 (for square lattice), which can be seen from figure 8 (cf. solid lines marked by triangles and crosses).

In the analysis above we have used the toy models in order to understand the spread of epidemics in systems of morphologically realistic roots. The models were designed in such a way as to capture the main features of root systems that influence the spread of epidemics. Our analysis of the toy models reveals that anisotropy, rather than disorder in root shape, has the major effect on epidemic invasion. The anisotropy in the real-root systems is due to the fact that the three-dimensional root density distribution is not centrally located beneath the seed. By mimicking a similar effect within toy models, we find that the greater the degree of root anisotropy, the more resistant the system is to epidemic invasion. The effect is due to anisotropy-induced, short-range correlations in disease transmissibilities between different root systems, which result in a reduced probability of invasion. These correlations are such that a root system that is well connected with one of its neighbours is likely to be poorly connected with its neighbour on the opposite side.

The other root characteristic, morphological complexity, has an opposite effect on the probability of epidemic invasion. By morphological complexity, we mean that the roots branch, producing secondary roots so that the shape of the root system becomes irregular. The increase in degree of branching makes the system more isotropic and diminishes the correlation effects, thus enhancing the vulnerability of the system to epidemics. This has been demonstrated with toy models by introducing disorder, characterized by the parameter *σ*_{J} mimicking the degree of branching. Overall, we can conclude that highly anisotropic and poorly branching root systems render the population less prevalent to invasion of soil-borne pathogens.

In our analysis, we have used the roots of young bean seedlings, grown under hydroponic conditions. These serve to illustrate the methodology. The hydroponically grown bean roots exhibit an inherent anisotropy in their shapes, which is the main morphological factor in reducing the invasion probability. Further work will consider the root morphology of monocotyledenous as well as dicotyledenous root systems grown and imaged in soil. We also assume quenched disorder in the systems considered here. This means that all the properties of the hosts, including their morphology, do not depend on time. Such an assumption is true only in the case of relatively fast epidemics and/or slowly growing root systems. Of course, in real situations, the root morphology can significantly change in the course of an epidemic, requiring time-dependent transmissibilities that are beyond the scope of this paper but constitute an important topic for further study.

## Acknowledgements

T.P.H. would like to thank the UK EPSRC for financial support. F.J.P.-R., S.N.T., F.M.N. and C.A.G. thank BBSRC for funding (grant no. BB/E017312/1). L.d.F.C. thanks CNPq (301303/2006-1 and 573583/2008-0) and FAPESP (05/00587-5) for sponsorship. M.M. thanks FAPESP (07/50988-1) for his grant and C.A.G. gratefully acknowledges the support of a BBSRC Professorial Fellowship. The authors thank Alexandre Cristino for suggestions regarding the growth of the seeds.

- Received June 2, 2010.
- Accepted July 5, 2010.

- © 2010 The Royal Society