Stochastic spreading models reproduce embolism propagation dynamics in angiosperm xylem networks of vessels connected by bordered pits

Onerva Korhonen; Steven Jansen; Luciano de Melo Silva; Petri Kiuru; Magdalena Held; Anna Lintunen; Annamari Laurén

doi:10.1017/nws.2026.10035

Stochastic spreading models reproduce embolism propagation dynamics in angiosperm xylem networks of vessels connected by bordered pits

Published online by Cambridge University Press: 25 June 2026

Onerva Korhonen

Steven Jansen ,

Luciano de Melo Silva ,

Anna Lintunen and

Onerva Korhonen*: Affiliation:
Data Science Research Centre, Tampere University , Tampere, Finland Tampere Institute for Advanced Study, Tampere University , Tampere, Finland School of Forest Sciences, University of Eastern Finland , Joensuu, Finland Department of Computer Science, Aalto University , Helsinki, Finland
Steven Jansen: Affiliation:
Institute of Botany, Ulm University, Ulm, Germany
Luciano de Melo Silva: Affiliation:
Institute of Biology, University of Graz, Graz, Austria
Petri Kiuru: Affiliation:
School of Forest Sciences, University of Eastern Finland , Joensuu, Finland
Magdalena Held: Affiliation:
Institute for Atmospheric and Earth System Research, University of Helsinki, Helsinki, Finland Department of Forest Sciences, Faculty of Agriculture and Forestry, University of Helsinki, Helsinki, Finland
Anna Lintunen: Affiliation:
Institute for Atmospheric and Earth System Research, University of Helsinki, Helsinki, Finland
Annamari Laurén: Affiliation:
School of Forest Sciences, University of Eastern Finland , Joensuu, Finland Department of Forest Sciences, Faculty of Agriculture and Forestry, University of Helsinki, Helsinki, Finland
*: Corresponding author: Onerva Korhonen; Email: onerva.korhonen@tuni.fi

Article contents

Abstract
Introduction
Materials and methods
Results
Discussion
Conclusion
Supplementary material
Competing interests
Funding statement
Data availability statement
Ethical standards
Author contributions
References

Rights & Permissions

Abstract

Plant xylem consists of a network of interconnected vessels, through which water is transported under negative pressure. Filling of vessels with air, or embolism, disturbs this transport process and, in extreme cases, leads to tree mortality. Despite this significance, embolism propagation dynamics are still poorly understood, primarily because xylem is opaque to direct observation. Furthermore, existing models of embolism spreading build excessively on physiological and anatomical parameters, and many misrepresent the intervessel pit membrane as a 2D surface. Here, we first extend these physiological models by implementing the pit membrane as a 3D object. Then, we introduce a susceptible-infected (SI) model, a simple stochastic model for tracking spreading through a population, for embolism propagation. After correctly fitting the spreading probability, our SI model reproduces vulnerability curves produced by both the physiological model and empirical data, highlighting that the SI model can address embolism spreading dynamics in plant species, for which detailed physiological data are not available. Furthermore, relating the SI model to the physiological one allows interpreting embolism spreading as a directed percolation process. Elucidating the exact mapping between directed percolation and embolism spreading will likely yield new fundamental insights into the relationships between xylem network architecture and embolism dynamics.

Keywords

Compartment models SI model epidemic spreading directed percolation plant physiology xylem hydraulics xylem networks embolism spreading Betula pendula Roth intervessel pit membrane

Information

Type: Research Article
Information: Network Science , Volume 14 , 2026 , e16

DOI: https://doi.org/10.1017/nws.2026.10035 [Opens in a new window]
Creative Commons: This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright: © The Author(s), 2026. Published by Cambridge University Press

1. Introduction

Water is transported in xylem tissue of plants as a continuous column through interconnected dead cells, known as vessels, under negative pressure. This transport is driven by the root-leaf pressure gradient caused by water evaporation at the leaves (Dixon and Joly, Reference Dixon and Joly1895). However, environmental conditions, in particular prolonged drought, can disturb this transport process by provoking embolism, or filling of xylem vessels with air. As embolism reduces the number of functional xylem vessels, sap transport capacity declines unless new vessels are developed (Tyree and Sperry, Reference Tyree and Sperry1989). Hydraulic failure in xylem, together with carbon starvation, is recognized as the major cause for tree mortality during drought (Adams et al., Reference Adams, Zeppel, Anderegg, Hartmann, Landhäusser, Tissue and Allen2017; Anderegg et al., Reference Anderegg, Klein, Bartlett, Sack, Pellegrini, Choat and Jansen2016; Choat et al., Reference Choat, Brodribb, Brodersen, Duursma, López and Medlyn2018; González-Muñoz et al., Reference González-Muñoz, Sterck, Torres-Ruiz, Petit, Cochard, von and Heinze2018; Mantova et al., Reference Mantova, Herbette, Cochard and Torres-Ruiz2022; McDowell et al., Reference McDowell, Pockman, Allen, Breshears, Cobb, Kolb and Williams2008; Sevanto et al., Reference Sevanto, Mcdowell, Dickman, Pangle and Pockman2014; Tavares et al., Reference Tavares, Oliveira, Mencuccini, Signori-Müller, Pereira, Diniz and Galbraith2023).

While the origin of vessel embolism is not fully known, one prominent mechanism is expansion of gas nanobubbles (Ingram et al., Reference Ingram, Reischl, Vesala and Vehkamäki2024) that appear in the sap more commonly than traditionally thought (Schenk et al., Reference Schenk, Steppe and Jansen2015). Embolism events are not randomly distributed among xylem vessels. Instead, embolism is known to spread from embolized vessels to their sap-filled neighbors (Avila et al., Reference Avila, Guan, Kane, Cardoso, Batz, DaMatta and McAdam2022; Choat et al., Reference Choat, Badel, Burlett, Delzon, Cochard and Jansen2016; Guan et al., Reference Guan, Pereira, McAdam, Cao and Jansen2021; Kaack et al., Reference Kaack, Weber, Isasa, Karimi, Li, Pereira and Jansen2021; Silva, Pereira et al., Reference Silva, Pereira, Kaack, Guan, Pfaff, Trabi and Jansen2024). However, the field of plant physiology still lacks a golden standard model for embolism spreading that would connect the microscopic function of pit membranes with the mesoscopic xylem structure.

One of the factors that affect embolism spreading in xylem is the interconnectivity between vessels (Loepfe et al., Reference Loepfe, Martinez-Vilalta, Piñol and Mencuccini2007). Xylem network models (e.g., Loepfe et al. (Reference Loepfe, Martinez-Vilalta, Piñol and Mencuccini2007); Mrad et al. (Reference Mrad, Domec, Huang, Lens and Katul2018); Mrad et al. (Reference Mrad, Johnson, Love and Domec2021)) address this interconnectivity by modeling the xylem as a complex network, or a system of nodes connected by links. Nodes of xylem network models represent the vessels, and links depict the bordered pits between neighboring vessels. These models build typically on meticulous information about xylem anatomy and physiology, including pit dimensions. For example, Mrad et al. (Reference Mrad, Domec, Huang, Lens and Katul2018) and Mrad et al. (Reference Mrad, Johnson, Love and Domec2021) introduced a model for constructing xylem networks with properties (e.g., average vessel length, vessel density) matching those observed empirically. While this model constructs the xylem network stochastically, the spreading of embolism between neighboring vessels is deterministic and depends on the bubble propagation pressure (BPP) dictated by intervessel pit membrane properties. Mrad et al. (Reference Mrad, Domec, Huang, Lens and Katul2018) and Mrad et al. (Reference Mrad, Johnson, Love and Domec2021) draw the BPPs from Weibull distribution defined by the largest 10% of measured pore diameters. However, this model does not account for the three-dimensionality of intervessel pit membranes that has important effects on embolism propagation (Kaack et al., Reference Kaack, Weber, Isasa, Karimi, Li, Pereira and Jansen2021; Pereira et al., Reference Pereira, Kaack, Guan, Silva, Miranda, Pires, Ribeiro, Schenk and Jansen2023; Zhang et al., Reference Zhang, Pereira, Kaack, Liu and Jansen2024). As the intervessel pit membrane pores are rather 3D tunnels through the membrane than 2D openings, the resistance of the pit membrane depends on the diameter of the smallest pore constriction instead of the diameter of the pore opening measured at the surface of the pit membrane. However, the 3D aspects of intervessel pit membranes as mesoporous media have not yet been implemented in vessel-level models. Here, we extend the Mrad model by implementing intervessel pits as 3D pit structures (Kaack et al., Reference Kaack, Weber, Isasa, Karimi, Li, Pereira and Jansen2021), where each pore pathway comprises several pore constrictions that together form a channel through the pit membrane.

Physiological models of embolism spreading (e.g., Kaack et al. (Reference Kaack, Weber, Isasa, Karimi, Li, Pereira and Jansen2021); Loepfe et al. (Reference Loepfe, Martinez-Vilalta, Piñol and Mencuccini2007); Lu et al. (Reference Lu, Smith-Martin, Muscarella, Uriarte and Zheng2025); Mrad et al. (Reference Mrad, Domec, Huang, Lens and Katul2018); Mrad et al. (Reference Mrad, Johnson, Love and Domec2021); Wason et al. (Reference Wason, Bouda, Lee, McElrone, Phillips, Shackel and Brodersen2021)) have made progress in explaining the microscale processes behind embolism spreading. However, the detailed anatomical and physiological data required for constructing these models are available only from a limited number of species, given the arduous measurements required for collecting those data. Further, the physiological models are sensitive to inaccurate assumptions about the complex multiphase interactions between gas, xylem sap with polar lipids, and intervessel pit membranes. To address this problem, we model embolism spreading as a stochastic process, taking place within a complex xylem network, in terms of the susceptible-infected (SI) spreading model. The SI model is a member of the compartment models family (Hethcote, Reference Hethcote2000; Pastor-Satorras et al., Reference Pastor-Satorras, Castellano, Van Mieghem and Vespignani2015). In network science, compartment models have been applied successfully to investigate diverse spreading phenomena, including infectious diseases (Barrat et al., Reference Barrat, Cattuto, Kivelä, Lehmann and Saramäki2021; Newman, Reference Newman2002; Rizi et al., Reference Rizi, Faqeeh, Badie-Modiri and Kivelä2022), rumor spreading (Nekovee et al., Reference Nekovee, Moreno, Bianconi and Marsili2007; Raponi et al., Reference Raponi, Khalifa, Oligeri and Di Pietro2022; Shah and Zaman, Reference Shah and Zaman2011), and innovation adoption (Duanmu and Chai, Reference Duanmu and Chai2025). In plant physiology, Roth-Nebelsick (Reference Roth-Nebelsick2019) and Torres-Ruiz et al. (Reference Torres-Ruiz, Cochard, Mencuccini, Delzon and Badel2016) have introduced compartment models for embolism spreading in a continuous medium. Compared to the physiological embolism spreading models, the main advantage of the SI model is its lightness: instead of the dozen physiological and anatomical parameters that characterize the physiological embolism spreading models, the SI model operates with a single parameter, namely spreading probability ( $\beta$ ).

Here, we investigate, to what extent the SI model can replicate embolism spreading produced by a physiological model or extracted from empirical data in the xylem network of Betula pendula Roth. We expect that after selecting an optimal spreading probability value for the SI model, the outcomes of the model are comparable to those produced by the physiological model or extracted from empirical data. We assume that this similarity applies both to the percentage loss of conductance (PLC) at the end of the spreading process and to the evolution of hydraulic and network-theoretical properties of the xylem network (e.g., the connected component size and fraction of embolized vessels) during the process. Further, we hypothesize that the vessel pathways in xylem break down faster than the total volume of embolized vessels increases. Thus, we expect that the hydraulic failure related to embolism spreading is primarily caused by xylem network breakdown, rather than the loss of sap in wood tissue. According to this hypothesis, embolism spreading can be interpreted as an instance of a universal network phenomenon, namely directed percolation.

2. Materials and methods

2.1 Xylem network model

We modeled the xylem anatomy with the network model introduced by Loepfe et al. (Reference Loepfe, Martinez-Vilalta, Piñol and Mencuccini2007), Mrad et al. (Reference Mrad, Domec, Huang, Lens and Katul2018), and Mrad et al. (Reference Mrad, Johnson, Love and Domec2021). This model represents the xylem as a 3D cylindrical grid of size $N_z \times N_r \times N_\phi$ with the periodic boundary condition $\phi _{N_\phi + 1} = \phi _0$ , that is, the number of independent angles considered was $N_\phi$ (Figure 1A–B). The row ( $z$ ), column ( $r$ ), and angular ( $\phi$ ) axes of the grid correspond to the axial, radial, and tangential directions of wood tissue. To create the xylem network, we considered the grid cells one by one in the order of their $\phi$ , $r$ , and $z$ coordinates. For each cell $(z, r, \phi )$ , we checked if the cell could continue an existing vessel, that is, if the corresponding cell at the previous row, $(z-1, r, \phi )$ belonged to a vessel. In this case, the cell $(z, r, \phi )$ ended the vessel (i.e., did not continue it) with probability $\mathcal{P}_c$ (for all abbreviations and symbols used in the present article, see Table 1). On the other hand, if $(z-1, r, \phi )$ did not belong to an existing vessel, the cell $(z, r, \phi )$ started a new vessel with probability $\mathcal{N}\mathcal{P}_c$ . There were two special cases: cells of the first row (i.e., $z = 0$ ) could not end a vessel and, similarly, cells of the last row (i.e., $z = N_z$ ) could not start a new vessel. Further, all vessels that reached the last row ended automatically independently of the $\mathcal{P}_c$ parameter. To define vessel diameters, we draw a set of values, corresponding to the number of vessels, from a lognormal distribution, the mean of which was equal to the empirically observed average vessel diameter of B. pendula (see Table 2). We assigned the diameters to vessels in relation to vessel length, that is, so that the longest vessel had the largest diameter. A Python implementation of our network construction model, together with the embolism spreading models and all analysis code used in the present article, is available at https://github.com/onerva-korhonen/tree-percolation.

Figure 1.

Schematic presentation of the xylem network and pore constriction models. A) The xylem is represented by a 3D cylindrical grid where the row ( $r$ ), column ( $z$ ), and angular ( $\phi$ ) dimensions correspond to the axial, radial, and tangential dimensions of a tree. Each grid cell can either belong to a vessel (red cells) or not (gray cells). Neighboring vessels form intervessel connections in radial direction with probability $\mathcal{P}_{e, rad}$ and in tangential direction with probability $\mathcal{P}_{e, tan}$ . B) In each angular slice of the xylem, each cell above a non-vessel cell starts a new vessel with probability $\mathcal{N}\mathcal{P}_c$ , while each cell above a vessel cell ends the existing vessel with probability $\mathcal{P}_c$ . Note that isolated vessels and dead-ends are removed before further analysis. C) Each pit membrane pore is modeled as a set of constrictions (red), forming a channel through the 3D pit membrane (gray). The smallest constriction diameter $R_{min}$ defines the effective diameter of each pore, and BPP of a pit membrane depends on the largest $R_{min}$ among all pores of the membrane.

A schematic representation of the xylem network and pore constriction models in plants.

Figure 1 Long description

Panel A: The xylem is depicted as a 3D cylindrical grid with axial, radial, and tangential dimensions. The grid cells are colored red to indicate vessels and gray to indicate no vessels. Neighboring vessels form intervessel connections in the radial and tangential directions with specified probabilities. Panel B: Each angular slice of the xylem shows cells above non-vessel cells starting new vessels with a certain probability, while cells above vessel cells end the existing vessel with another probability. Isolated vessels and dead-ends are removed before further analysis. Panel C: Each pit membrane pore is modeled as a set of constrictions forming a channel through the 3D pit membrane. The smallest constriction diameter defines the effective diameter of each pore, and the largest pore diameter determines the bubble point pressure of the pit membrane.

This procedure yielded a set of vessels with known physical locations. These vessels formed the nodes of the xylem network, while the links of the network corresponded to the connections between neighboring vessels, known as vessel multiples (VMs). To create these links, we identified the neighboring cells, or cells sharing a face or a vertex in $r$ or $\phi$ directions, of each cell belonging to a vessel. A VM connected two neighboring cells in different vessels with probability $\mathcal{P}_{e, rad}$ in the radial direction and with probability $\mathcal{P}_{e, tan}$ in the tangential direction.

Several VMs can exist between a vessel pair. We considered all VMs when calculating the sap flow through the xylem and the related effective conductance of the xylem network (Section 2.2). However, for the sake of simplicity, in the SI model for embolism spreading (Section 2.3.2), we considered the vessel network unweighted, meaning that all links have an equal weight, independently of the pit area shared by the vessels connected by the link.

The procedure for creating vessels and VMs may lead to isolated vessels. Such vessels are rare or nonexistent in real plants (Bosshard and Kučera, Reference Bosshard and Kucera1973; Brodersen et al., Reference Brodersen, Lee, Choat, Jansen, Phillips, Shackel and Matthews2011; Burggraaf, Reference Burggraaf1972). Thus, we considered them artifacts of our network construction algorithm and removed them from any further analysis. Similarly, we removed any dead-ends, that is, vessel groups that were not connected to both the inlet and the outlet rows, following the example of earlier xylem network construction algorithms (Brodersen et al., Reference Brodersen, Lee, Choat, Jansen, Phillips, Shackel and Matthews2011; Mrad et al., Reference Mrad, Domec, Huang, Lens and Katul2018). Importantly, our model included only one type of conducting cells, vessels.

Table 1.

Overview of the abbreviations and symbols used. For values of parameters, see Table 2

A table listing abbreviations, units, and definitions related to xylem network modeling.

Table 2.

Xylem network and embolism spreading model parameters for B. pendula

A table listing parameters for a xylem network and embolism spreading model.

^aVessel diameter values are drawn from a lognormal distribution with the given mean and assigned to vessels so that the longest vessel gets the largest diameter.

^bObtained by multiplying the vessel density from Lintunen and Kalliokoski (Reference Lintunen and Kalliokoski2010), 630.82 1/mm², by the average vessel crosssection area.

^cEstimated based on Table 1 of Kaack et al. (Reference Kaack, Weber, Isasa, Karimi, Li, Pereira and Jansen2021) using $T_m$ and $t_f$ and assuming 10 nm between microfibril strands.

2.2 Hydraulic properties of the xylem network

We measured the functionality of the xylem network in terms of effective conductance $K_{eff}$ defined as the total sap flow through the network normalized by the pressure difference between the network inlet and outlet. Following Mrad et al. (Reference Mrad, Domec, Huang, Lens and Katul2018), we calculated the sap flow $Q_I$ through vessel $I$ in terms of the Hagen–Poiseuille equation:

(1)

\begin{equation} Q_I = \frac {\pi D_I^4}{128\mu L_{v,I}}\Delta P_I, \end{equation}

where $D_I$ denotes the diameter and $L_{v,I}$ the length of vessel $I$ , $\mu$ is sap viscosity, and $\Delta P_I$ is the pressure difference between the ends of the vessel. Further, the flow through VM $(I,J)$ is obtained by considering the Sampson flow resistance for pores in a plate and a Poiseuille flow contribution through the finite membrane thickness as:

(2)

\begin{equation} Q_{(I,J)} = \frac {D_p^3N_{p,(I,J)}}{24\mu }\left (1 - \frac {16T_m}{3\pi D_p} \right )^{-1}\Delta P_{(I,J)}, \end{equation}

where $D_p$ is the average intervessel pit membrane pore diameter, $T_m$ is the average intervessel pit membrane thickness, and $\Delta P_{(I,J)}$ is the pressure difference over the bordered intervessel pit pair (Dagan et al., Reference Dagan, Weinbaum and Pfeffer1982; Mrad et al., Reference Mrad, Domec, Huang, Lens and Katul2018). The number of pores $N_{p,(I,J)}$ is obtained, following Mrad et al. (Reference Mrad, Domec, Huang, Lens and Katul2018) and Sperry and Hacke (Reference Sperry and Hacke2004), by dividing the VM pit field area $A_{I,J}$ by the area of an average pit membrane pore:

(3)

\begin{equation} N_{p,(I,J)} = \frac {A_{(I,J)}}{\pi \left (D_p + t_f \right )^2/4}, \end{equation}

where $t_f$ is microfibril strand thickness. The formula for $A_{(I,J)}$ is given by Mrad et al. (Reference Mrad, Domec, Huang, Lens and Katul2018) and Loepfe et al. (Reference Loepfe, Martinez-Vilalta, Piñol and Mencuccini2007) as

(4)

\begin{equation} A_{(I,J)} = \frac {1}{2}\left (\frac {A_{I}}{N_{VM,I}} + \frac {A_J}{N_{VM,J}}\right )f_cf_{pf}. \end{equation}

Here, the fraction of the surface area $A_I$ and number of VMs $N_{VM,I}$ of vessel $I$ defines the theoretical vessel surface area available for VM $(I,J)$ , the vessel contact fraction $f_c$ tells the average fraction of wall area shared with other vessels out of total vessel surface area, and $f_{pf}$ is the fraction of the contact area between two vessels, $A_If_c$ , occupied by intervessel pit membranes.

Combining equations 1 and 2 allowed obtaining the total water flow $Q$ through the system, assuming continuity everywhere except for the inlet and outlet and given the inlet and outlet pressures $P_{in}$ and $P_{out}$ . Finally, from $Q$ we obtained the effective conductance as

(5)

\begin{equation} K_{eff} = \frac {Q}{P_{in} - P_{out}}. \end{equation}

As stated above, the total flow $Q$ , and consequently $K_{eff}$ calculated through equation 5, originate from a combination of Pouseuille and Sampson flows. Both these flows are solutions to the general Navier–Stokes equations for the motion of viscous fluids derived under special conditions. Therefore, in our Python implementation, we simulated the sap flow using the StokesFlow algorithm of the OpenPNM Python package (Gostick et al., Reference Gostick, Aghighi, Hinebaugh, Tranter, Hoeh, Day and Putz2016), version 3.4.1.

2.3 Embolism spreading

To investigate the vulnerability of the xylem network to embolism, we simulated embolism spreading between vessels with two approaches. First, we considered a physiological spreading model based on the pit-level properties of VMs. Second, we showed that the SI model (Hethcote, Reference Hethcote2000; Pastor-Satorras et al., Reference Pastor-Satorras, Castellano, Van Mieghem and Vespignani2015) can replicate the results of the physiological spreading model while significantly reducing the amount of free model parameters and thus reducing the model complexity. In both models, as soon as a vessel became embolized, it was removed from the network together with all its VMs, so that it no longer contributed to sap flow or any network properties (see Section 2.4). In other words, at any given moment, the intervessel network contained only sap-filled vessels, and thus the number of nodes in the network decreased during the embolism spreading process.

Both models assumed that there is at least one embolized vessel, the spreading seed, in the xylem at the beginning of the embolism spreading process. This assumption is reasonable for xylem tissue of perennial, woody plants because embolism is often observed in older growth rings or primary xylem (Choat et al., Reference Choat, Badel, Burlett, Delzon, Cochard and Jansen2016). Our models described how this initial embolism spreads in the xylem but did not aim to explain mechanistically how the seed vessel became embolized. Further, we did not consider refilling of vessels as the time scale of refilling is notably slower than the embolism spreading time scales we considered here. Further, while seasonal refilling takes place in B. pendula under positive pressures (Strati et al., Reference Strati, Patiño, Slidders, Cundall and Mencuccini2003), this recovery mechanism is less probable under normal, negative xylem pressures (Choat et al., Reference Choat, Brodribb, Brodersen, Duursma, López and Medlyn2018).

2.3.1 Physiological spreading model

In the physiological spreading model, the spreading of embolism between two vessels depended on the VM-specific BPP. We calculated the BPPs using the pore constriction model introduced by Kaack et al. (Reference Kaack, Weber, Isasa, Karimi, Li, Pereira and Jansen2021) (referred to as Model 1 in the original article) (Figure 1C). In other words, we modeled the intervessel pit membrane as a 3D object so that each pit membrane pore formed a channel, consisting of several constrictions, through the 3D pit membrane. We drew the constriction diameters from a left-truncated normal distribution with mean value $\bar {r}$ , standard deviation $\sigma _r$ , and cutoff $R_L$ , and the smallest constriction diameter defined the diameter of the pore. The number of constrictions per pore depended on the intervessel pit membrane thickness as

(6)

\begin{equation} N_{pc} = (T_m + 20) / 30, \end{equation}

where $T_m$ was expressed in nanometers and a 10 nm spacing between microfibril strands was assumed (Kaack et al., 2021). To define the BPP of an intervessel pit field with area $A_p$ , we simulated a set of pores so that their total area was as close as possible to $A_p$ without exceeding it. Calculation of the pit BPP was based on the Young–Laplace equation, a general equation that describes how the pressure difference across the interface between two fluids depends on the interface curvature and surface tension. Following Kaack et al. (Reference Kaack, Weber, Isasa, Karimi, Li, Pereira and Jansen2021) and Schenk et al. (Reference Schenk, Steppe and Jansen2015), we calculated BPP using a modified Young–Laplace equation as

(7)

\begin{equation} \textrm {BPP} = \kappa 2 \gamma cos(\alpha )/R_{min, max}, \end{equation}

where $\kappa$ is a dimensionless pore shape correction factor that accounts for irregularities in pore shape, in particular for deviation from the ideal cylindrical shape, for which the modified Young–Laplace equation has originally been derived (Emory, Reference Emory1989; Schenk et al., Reference Schenk, Steppe and Jansen2015; US EPA Office of Water, 2005). Further, $\alpha$ is the gas contact angle between the propagating air bubble and the sap surface and $R_{min, max}$ the largest pore diameter, that is, the largest among the smallest constriction diameters found in all the pores of the intervessel pit membrane. Finally, $\gamma$ is the equilibrium value of the dynamic sap surface tension that depends on the local lipid concentration at sap–air interfaces; this value is notably smaller than the surface tension of pure water (Kaack et al., Reference Kaack, Altaner, Carmesin, Diaz, Holler, Kranz and Donaldson2019) (for exact parameter values used, see Table 2). For each pair of neighboring vessels $I$ and $J$ , we defined BPP( $I,J$ ) as the smallest of the BPPs of the VMs connecting $I$ and $J$ .

When defining the number of pores in a pit membrane, we calculated the area of each pore using the smallest constriction diameter as pore diameter. This gives a theoretical upper limit for the number of pores for a case where all constrictions of each pore have the same, minimum diameter. Similarly, calculating pore area based on the largest constriction diameter gives a theoretical lower limit that corresponds to a case where the largest constrictions of each pore are aligned. In reality, pores have constrictions larger than the smallest one, constrictions of different pores are not aligned, and the pore tunnels may be curved, which makes estimating the actual number of pores per pit challenging. While we assume that the effects of estimating the number of pores in terms of the theoretical upper limit to be small, this decision may have yielded slightly underestimated BPPs and, correspondingly, lower embolism resistance and less negative $P_{50}$ values.

Air moves from an embolized vessel to its sap-filled neighbors through a series of snap-off events creating nanobubbles. These bubbles may expand and embolize the sap-filled neighboring vessel. On the other hand, they may become stabilized, in which case their dynamic surface tension, caused by the stretching of the surfactant molecules coating the bubble, prevents expansion (Kaack et al., Reference Kaack, Altaner, Carmesin, Diaz, Holler, Kranz and Donaldson2019; Schenk et al., Reference Schenk, Steppe and Jansen2015; Silva et al., Reference Silva, Bujnowski, Pereira, Miranda, Schenk and Jansen2025). The fracture of this coating typically leads to sudden increase of surface tension and collapse of the nanobubble (Kaack et al., Reference Kaack, Altaner, Carmesin, Diaz, Holler, Kranz and Donaldson2019; Schenk et al., Reference Schenk, Steppe and Jansen2015; Silva et al., Reference Silva, Bujnowski, Pereira, Miranda, Schenk and Jansen2025). However, modeling the coating process would require detailed information about the multiphase interactions inside the intervessel pit membrane, including the nanoscale interactions between surfactants, xylem sap, and gas bubbles. For simplicity, we did not model this process explicitly but instead assumed that embolism spreading happens between vessels $I$ and $J$ every time when the sap pressure exceeds the BPP. Therefore, at each time step $t$ of the spreading simulation, we embolized each non-embolized vessel $I$ that had an embolized neighbor $J$ (i.e., a neighbor $J$ with embolization time $\lt t$ ) so that BPP( $I, J$ ) $\le P$ , where $P$ was the sap pressure in the xylem.

2.3.2 SI model

The SI model is a stochastic model that enables tracking spreading through a population straightforwardly and with a minimum number of free parameters. This model and related compartment models are widely applied to study diverse phenomena including spreading of pathogens (e.g., Barrat et al. (Reference Barrat, Cattuto, Kivelä, Lehmann and Saramäki2021); Rizi et al. (Reference Rizi, Faqeeh, Badie-Modiri and Kivelä2022)), propagation of rumors and fake news (e.g., Nekovee et al. (Reference Nekovee, Moreno, Bianconi and Marsili2007); Shah and Zaman (Reference Shah and Zaman2011)), innovation adoption (e.g., Duanmu and Chai (Reference Duanmu and Chai2025)), and spreading of invasive species (e.g., Strickland et al. (Reference Strickland, Dangelmayr, Shipman, Kumar and Stohlgren2015)). The SI model divides the population (here, the xylem vessels) into two compartments: the susceptible (S) compartment, which is vulnerable to an infection but not yet infected, and the infected (I) compartment that contains infected individuals. At each time step, every individual in the S compartment can catch the infection from any of their neighbors in the I compartment with probability $\beta$ . Since these spreading events are independent, the total infection probability of a susceptible individual with $N_{inf}$ infected neighbors is

(8)

\begin{equation} P_{inf} = 1 - (1 - \beta )^{N_{inf}}. \end{equation}

This S-I transition is the only transition between compartments in the SI model. The model does not contain any recovery mechanisms, and thus transitions from the I compartment to the S compartment are not possible.

In our SI model of embolism spreading, the S compartment represented the sap-filled vessels, while the I compartment corresponded to the embolized ones.

2.3.3 Vulnerability curves

To investigate the xylem network’s resistance to embolism under different conditions, we first simulated 100 individual xylem networks starting from a 100 × 10 × 100 grid (i.e., each network has at most 100 × 10 × 100 = 100,000 vessel elements). For each of these networks, we selected the seed vessels, from which the embolism spreading started, by identifying the components (i.e., groups of nodes that were reachable from each other through the network links) of the xylem network and picking a random vessel per component. These seed vessels represent the embolized vessels that are present in real-world xylem (Carmesin et al., Reference Carmesin, Port, Böhringer, Gottschalk, Rasche and Jansen2023).

For each of these networks, we simulated embolism spreading with the physiological spreading model for a set of 79 $P$ values: 8 values ranging from 0 MPa to −0.875 MPa with the step size of 0.125 MPa, 50 values ranging from −1 MPa to −1.49 MPa with the step size of 0.01 MPa, and final 21 values ranging from −1.5 MPa to −4 MPa with the step size of 0.125 MPa. We selected the uneven step size to emphasize the pressure area around the expected $P_{50}$ value (i.e., the pressure where 50% of $K_{eff}$ is lost). We set the embolization time of the seed vessels to zero and ran the physiological spreading model for each $P$ until the spreading saturated, that is, until the fraction of embolized vessels did not change for 20 time steps. The final $K_{eff, phys}$ for each $P$ was averaged across the 100 simulation iterations.

Then, for each $P$ , we calculated the PLC by comparing the $K_{eff, phys}$ at spreading saturation point to $K_{eff, 0}$ calculated for the xylem with no embolized vessels. Finally, we constructed the physiological vulnerability curve by drawing the PLC as a function of $P$ .

To find the SI model spreading probability $\beta$ that best approximates embolism spreading at each $P$ , we ran for each of the simulated networks the SI model for 65 $\beta$ values: 30 values ranging from 0 to 0.15 with step size of 0.005 and additional 35 values ranging from 0.175 to 1 with step size of 0.025. We selected the uneven step size through preliminary simulations that used an evenly spaced probability range from 0 to 1 with step size of 0.15. The outcomes of these simulations hinted that lower $\beta$ values yield a better fit to the physiological spreading model. Thus, we selected for the final analysis an uneven step size that emphasizes the lower $\beta$ values.

Again, we ran the model until the spreading saturated and calculated the final $K_{eff, SI}$ as an average across the 100 iterations. Then, we selected for each $P$ the $\beta$ that produced $K_{eff, SI}$ closest to the corresponding $K_{eff, phys}$ . Next, we calculated the PLC for the selected $\beta$ values and drew the SI-based vulnerability curve as a function of $P$ .

In addition to optimizing $\beta$ against the physiological model, we fitted the SI model also against an empirical B. pendula vulnerability curve published by González-Muñoz et al. (Reference González-Muñoz, Sterck, Torres-Ruiz, Petit, Cochard, von and Heinze2018). Unlike for the physiological model, the absolute $K_{eff}$ values were not available for the empirical data. Therefore, we optimized $\beta$ directly in terms of PLC, selecting for each $P$ value, separately, the $\beta$ value that produced the PLC closest to that extracted from the empirical data.

2.4 Network properties

We explored the embolism spreading in xylem further in terms of a set of network properties calculated at each time step of the spreading process. For all these properties, we reported the mean and standard deviation across 100 iterations. Note that as described in Section 2.3, we removed vessels from the intervessel network as soon as they became embolized, and thus only sap-filled vessels contributed to the network properties.

Largest connected component (LCC) size. Depending on the values of $\mathcal{P}_c$ , $\mathcal{N}\mathcal{P}_c$ , $\mathcal{P}_{e, rad}$ , and $\mathcal{P}_{e, tan}$ , the xylem network could be divided into several components. Every vessel of a component was reachable from other vessels of the same component following network links, while there were no links between vessels belonging to different components. The LCC size was defined as the number of vessels in the largest network component. Because of the systematical removal of embolized vessels from the network, the LCC contained sap-filled vessels only.

Functional LCC size. We considered a network component functional if it contained at least one vessel starting from the inlet row and at least one vessel ending at the outlet row. In other words, a functional component spanned the whole xylem space in vertical direction and thus contributed to the sap flow through the xylem and $K_{eff}$ . We defined functional LCC size as the number of vessels in the largest functional component. The difference between the LCC and the functional LCC was the sap flow: while both the LCC and the functional LCC consisted of sap-filled vessels, the LCC did not necessarily connect to both network inlet and outlet and thus there was no guaranteed sap flow through it.

Prevalence is an epidemiological measure defined as the fraction of infected individuals out of the total population. Here, we calculated prevalence as the number of embolized vessels divided by the total number of vessels in the intact xylem network.

Nonfunctional sap volume. Unlike the functional components, nonfunctional xylem network components lacked either the inlet, the outlet, or both. Thus, there was no sap flow through the nonfunctional components. These components arose from embolism spreading: while all components of a xylem network with no embolized vessels were functional and filled with sap, a component could lose its connection to the inlet or the outlet because of embolism. In such a case, the sap in the now nonfunctional component could not exit the xylem because of the lack of either access to outlet or the driving pressure difference between the ends of the component. We defined nonfunctional sap volume as the total volume of this sap trapped in the xylem, which was equal to the total volume of vessels in nonfunctional components.

Number of inlets and outlets per component was defined as the average number of vessels starting from the inlet row or ending at the outlet row per component.

The investigated network measures were related to each other. In particular, assuming that embolism spread at the same speed in all network components, we expected the relationship between the LCC size and the functional LCC size to be positive and approximately linear so that the LCC size set the upper boundary for the functional LCC size. We expected both the LCC size and the functional LCC size to decrease during embolism spreading, both due to the removal of embolized vessels from the network and because this removal could divide existing network components into several smaller ones. Further, we expected prevalence to increase approximately linearly with the decreasing LCC size, assuming that embolism events were equally distributed across network components. The expected relationship of the nonfunctional sap volume with the largest LCC size was more complicated. The nonfunctional sap volume increased with the decreasing LCC size when more components lost their connection to the inlet or to the outlet. However, the nonfunctional sap volume could also decrease during embolism spreading, as also vessels that belong to nonfunctional components could become embolized. Finally, we expected the number of inlets and outlets per component to decrease with the decreasing LCC size because of both the embolization of inlet and outlet vessels and the increased number of components. Importantly, all relationships between network measures depended on the intervessel connectivity structure and, for example, differences in connectivity structure between network components could lead to nonlinear relationships.

We evaluated the similarity of the behavior of network properties between the physiological and SI spreading models in terms of the root mean square error (RMSE) between the time-dependent behavior vectors. RMSE is defined for vectors of equal length. However, due to the different time scales of the two models, the SI model produced longer behavior vectors than the physiological model (see Section 3.2 for details). We addressed this by interpolating more values for the physiological model vector with the numpy.interp function and normalized the obtained RMSEs by the mean of the physiological behavior vector for easier comparison.

2.5 Parameter estimation

The xylem network model of Mrad et al. (Reference Mrad, Domec, Huang, Lens and Katul2018) and Mrad et al. (Reference Mrad, Johnson, Love and Domec2021) was parametrized for Acer glabrum. However, we selected B. pendula for the example species instead because this species shows a very broad geographic distribution pattern across the Eurasian continent. To re-parametrize the model for B. pendula, we used a set of parameters reported in the literature for large B. pendula branches (see Table 2).

To estimate the values of the probabilities $\mathcal{N}\mathcal{P}_c$ , $\mathcal{P}_c$ , $\mathcal{P}_{e, rad}$ , and $\mathcal{P}_{e, tan}$ , we constructed a set of xylem networks with the same size as those used for constructing the vulnerability curves so that there were 10 networks constructed with each combination of the probability parameters. We calculated for each network the vessel tissue fraction (VTF) (i.e., the fraction of the xylem cross-sectional area occupied by vessels) and vessel grouping index (GI) (i.e., the total number of vessels divided by the number of groups of adjacent vessels) and averaged them across the 10 networks per probability combination. We then compared the observed VTF and GI values to the values reported for B. pendula (Figure 2, see also Table 2). To this end, we assigned each parameter combination with two ranks, corresponding to VTF and GI, based on the absolute-valued difference between the value produced by the combination and the value reported in the literature. In other words, the parameter combination that produced the VTF value closest to the reference value from the literature, $\mathrm{VTF}_{ref}$ , received rank 1, the combination that produced the second-closest value received rank 2, and so on. We then added together the VTF and GI ranks and selected the parameter combination corresponding to the lowest total rank. This optimization yielded $\mathcal{N}\mathcal{P}_c = 0.8$ , $\mathcal{P}_c = 0.9$ , $\mathcal{P}_{e, tan} = 0.9$ , and $\mathcal{P}_{e, rad} = 0.0$ . However, we did not want to completely exclude the possibility of radial VMs. Thus, we used $\mathcal{P}_{e, rad} = 0.01$ instead of the optimization outcome. This parameter combination yielded VTF 0.26 and GI 2.61 (average across the 10 networks created for the selected probability parameter combination).

Figure 2.

Optimization of probability parameters for xylem network construction. We constructed xylem vessel networks with different combinations of probabilities $\mathcal{N}\mathcal{P}_c$ , $\mathcal{P}_c$ , $\mathcal{P}_{e, rad}$ , and $\mathcal{P}_{e, tan}$ and calculated the vessel tissue fraction (A, B) and grouping index (GI) (C, D) of each simulated network. Two 2D projections of the 4D space spanned by the four probability parameters are shown: for visualizing the effect of $\mathcal{N}\mathcal{P}_c$ and $\mathcal{P}_c$ (A, C), $\mathcal{P}_{e, rad}$ and $\mathcal{P}_{e, tan}$ are fixed to their optimal value, while the effect of $\mathcal{P}_{e, rad}$ and $\mathcal{P}_{e, tan}$ (B, D) is visualized fixing $\mathcal{N}\mathcal{P}_c$ and $\mathcal{P}_c$ to their optimal values. The red square shows the optimal parameter combinations. Note that the z scale of GI (C, D) is logarithmic; the blank cells indicate areas where GI = 0.

Figure 2 Long description

Four heatmaps depict the optimization of probability parameters for xylem network construction. Panel A: A heatmap shows the vessel tissue fraction as a function of P_c and N_Pc. The x-axis is labeled P_c, and the y-axis is labeled N_Pc. The color scale on the right indicates the vessel tissue fraction, ranging from 0.0 to 1.0. A red square marks the optimal parameter combination. Panel B: A heatmap displays the vessel tissue fraction as a function of P_e, tan and P_e, rad. The x-axis is labeled P_e, tan, and the y-axis is labeled P_e, rad. The color scale on the right indicates the vessel tissue fraction, ranging from 0.0 to 1.0. A red square marks the optimal parameter combination. Panel C: A heatmap illustrates the grouping index (GI) as a function of P_c and N_Pc. The x-axis is labeled P_c, and the y-axis is labeled N_Pc. The color scale on the right indicates the grouping index, ranging from 10⁰ to 10². A red square marks the optimal parameter combination. Panel D: A heatmap shows the grouping index (GI) as a function of P_e, tan and P_e, rad. The x-axis is labeled P_e, tan, and the y-axis is labeled P_e, rad. The color scale on the right indicates the grouping index, ranging from 10⁰ to 10². A red square marks the optimal parameter combination. The z scale of GI is logarithmic, and blank cells indicate areas where GI equals 0.

Network science commonly quantifies network structure in terms of average degree, defined as the average number of neighbors per network node, and network density, defined as

(9)

\begin{equation} D = \frac {2E}{N(N-1)}, \end{equation}

where $E$ is the number of network links (here, number of VMs) and $N$ is the number of nodes (here, number of vessels). These measures offer a basis for comparing the structure of networks of different sizes. The optimized probability parameters produced average degree 2.92 $\pm$ 2.87E-2 and network density 8.56E-4 $\pm$ 2.24E-5 (mean $\pm$ across 100 networks constructed with the optimized probability parameters and used for simulating embolism spreading). This indicated that the simulated intervessel networks were sparse, similar to many real-world networks (e.g., Watts and Strogatz (Reference Watts and Strogatz1998) reported $D = 2.71\mathrm{E-}4$ for a collaboration network of film actors, $D = 5.41\mathrm{E-}4$ for the power grid of the western US but notably more dense $D = 4.98\mathrm{E-}2$ for the neuronal network of Caenorhabditis elegans). Further, the low standard deviation of the degree suggested that the VMs were relatively equally distributed across the vessels.

3. Results

3.1 SI model replicated physiological embolism spreading

To investigate if the stochastic SI model can approximate spreading of embolism in B. pendula xylem, we simulated embolism spreading with two models. Our physiological model was based on xylem anatomical and physiological characteristics as well as the network models of Loepfe et al. (Reference Loepfe, Martinez-Vilalta, Piñol and Mencuccini2007), Mrad et al. (Reference Mrad, Domec, Huang, Lens and Katul2018), and Mrad et al. (Reference Mrad, Johnson, Love and Domec2021). The SI model, on the other hand, is a stochastic spreading model used to study epidemic spreading computationally in various fields, including network science (Hethcote, Reference Hethcote2000; Pastor-Satorras et al., Reference Pastor-Satorras, Castellano, Van Mieghem and Vespignani2015). To compare the two models, we selected the free parameter of the SI model, the spreading probability $\beta$ , separately for each $P$ to minimize the difference in effective conductance $K_{eff}$ between the two models. After this fitting, the SI model produced a vulnerability curve notably similar to the one produced by the physiological model (Figure 3A), demonstrating that the spreading saturated at similar levels of embolism in both models. The optimized $\beta$ value varied between 0 and 0.085, increasing almost monotonically as a function of $P$ (Figure 3B).

Figure 3.

SI model replicates the vulnerability curve produced by the physiological model. A) PLC produced by the physiological (red) and SI (gray) models. For calculating the PLC, $K_{eff}$ was averaged across 100 spreading simulations. The transition from normal xylem function to hydraulic failure happened almost step-like over a narrow $P$ range; for further details, see the main text and Section 4.3. B) The optimized SI spreading probability ( $\beta$ ).

Two line graphs compare xylem pressure effects on PLC simulated with the physiological and SI spreading models and the optimized beta.

The shape of the vulnerability curves was sigmoidal, almost step-like, indicating a sharp transition between normal xylem function and hydraulic failure. The $P_{12}$ , $P_{50}$ , and $P_{88}$ values (i.e., pressures where 12%, 50%, and 88% of $K_{eff}$ were lost) of our physiological model were −1.26 MPa, −1.28 MPa, and −1.29 MPa, respectively. The full spreading process for $P_{12}$ , $P_{50}$ , and $P_{88}$ is shown in Figures 4 and 5.

Figure 4.

Behavior of xylem network properties during embolism spreading simulated with the physiological model. $K_{eff}$ and the largest connected component (LCC) size (A, B, C), prevalence (D, E, F), and nonfunctional sap volume (G, H, I) are monitored for $P_{12}$ (A, D, G), $P_{50}$ (B, E, H), and $P_{88}$ (C, F, I) of the physiological spreading model. For each network property, the full line corresponds to the mean across 100 simulations, while the shadowed area shows the standard deviation across iterations. The x axis shows the dimensionless number of simulation time steps that does not directly map to absolute time. Note that the embolism spreading took different number of time steps to saturate at different $P$ values, and the x axes thus differ between subplots.

Multiple graphs depict the behavior of xylem network properties during embolism spreading simulated with a physiological model.

Figure 4 Long description

The figure contains nine line graphs organized in a 3x3 grid, labeled A to I. Each column represents different conditions with varying pressure values: P12, P50, and P88. Panel A, B, and C show the effective hydraulic conductivity (K_eff) and the largest connected component (LCC) size over time. The y-axis on the left measures K_eff in units of 1E-16 x m² x Pa^-1 x s^-1, and the y-axis on the right measures LCC size. The x-axis represents time in dimensionless simulation time steps. Panel D, E, and F display the prevalence of embolism over time, with the y-axis measuring prevalence and the x-axis representing time. Panel G, H, and I illustrate the nonfunctional sap volume over time, with the y-axis measuring volume in units of 1E-9 x m³ and the x-axis representing time. Each graph shows the mean across 100 simulations as a solid line, with the shadowed area indicating the standard deviation across iterations. The x-axes differ between subplots as the embolism spreading took different numbers of time steps to saturate at different values.

3.2 Network properties revealed the embolization process

Vulnerability curves showed the final state of embolism spreading at each $P$ . To explore how this final state was reached, we monitored several properties of the xylem network during the spreading process. In particular, we explored the behavior of network properties as a function of increasing number of simulation time steps. At each time step, an embolized vessel could spread the embolism to all its sap-filled neighbors. In other words, the time step length corresponded to the time required by a sap-filled vessel to become fully embolized and able to spread the embolism further. This time depends on both xylem properties and external conditions, and therefore mapping simulation time steps to absolute time is not trivial.

In the case of the physiological spreading model (Figures 4, S1A–B), the behavior of $K_{eff}$ , the LCC size, and the functional LCC size at $P_{50}$ and $P_{88}$ was characterized by a quick decrease followed by a long, slower saturation phase, where the standard deviation across simulations increased. On the other hand, at $P_{12}$ , $K_{eff}$ , the LCC size, and the functional LCC size stayed stable over time. In the beginning of the simulation, the LCC was functional for all $P$ s. However, for $P_{50}$ and $P_{88}$ , the size of the largest nonfunctional component exceeded the size of the largest functional component in the beginning of the saturation phase (Figure S1). Prevalence, or the fraction of embolized vessels, reflected the behavior of $K_{eff}$ : quick initial increase was followed by a long saturation phase at $P_{50}$ and $P_{88}$ , while the behavior at $P_{12}$ was approximately static. Interestingly, at $P_{88}$ , the average final prevalence was lower than the final PLC, while at $P_{12}$ and $P_{50}$ , the final prevalence and PLC were approximately equal. Behavior of nonfunctional sap volume, defined as the total amount of water in vessels with no connection to either inlet, outlet, or both, supported the observation about the nonfunctionality of the LCC at more negative pressures: at $P_{50}$ and $P_{88}$ , the nonfunctional sap volume increased as the embolism spreading proceeded and broke down network components. Finally, the average number of inlets and outlets showed the familiar pattern of quick initial decrease and slower saturation at $P_{50}$ and $P_{88}$ and static behavior at $P_{12}$ (Figure S2). The number of outlets was higher than that of inlets because of the xylem network construction process: while a cell in the inlet row initiated a vessel with probability $\mathcal{N}\mathcal{P}_c$ , all vessels that reached the outlet row ended automatically, thus yielding a higher number of outlet than inlet cells.

The most notable difference in the behavior of network properties between the physiological model and the SI model (Figures 5 and S1C–D) was the time scale. While the physiological spreading model saturated within less than 100 steps for $P_{12}$ , around 250 steps for $P_{50}$ , and around 140 steps for $P_{88}$ , the SI model necessitated over 3000 steps for $P_{12}$ , over 1000 steps for $P_{50}$ , and around 600 steps for $P_{88}$ . Despite the different time scales and the related rather high normalized RMSE values (Table 3), the behavior of $K_{eff}$ , the LCC size and the functional LCC size, prevalence, and the number of inlets and outlets per component were qualitatively similar between the two models, although the SI model yielded higher standard deviation at $P_{12}$ . The nonfunctional sap volume, instead, behaved differently in the SI model: it increased only in the later states of the simulation for both $P_{50}$ and $P_{88}$ , and for $P_{50}$ , its final value remained lower than in the physiological model. However, for $P_{88}$ , the nonfunctional sap volume behavior of the physiological model approached that of the SI model.

Table 3.

RMSE between the network property behavior vectors calculated with the physiological and SI spreading models

A table comparing RMSE values for different network properties between physiological and SI spreading models.

^aRMSE and normalized RMSE; RMSE was normalized by the mean of the physiological behavior vector.

Figure 5.

Behavior of xylem network properties during embolism spreading simulated with the SI model. Similarly as in Figure 4, $K_{eff}$ and the largest connected component (LCC) size (A, B, C), prevalence (D, E, F), and nonfunctional sap volume (G, H, I) are monitored for $P_{12}$ (A, D, G), $P_{50}$ (B, E, H), and $P_{88}$ (C, F, I), the full lines and shadowed areas show mean and standard deviation across simulations, and the x axis shows the number of simulation time steps. Similarly as in Figure 4, note the difference in the x axis maximum values and spacing between subplots.

Multiple graphs depict the behavior of xylem network properties during embolism spreading simulated with the SI model.

Figure 5 Long description

Here, subgraphs A, D and G correspond to P_12, subgraphs B, E, H correspond to P_50 and subgraphs C, F, I correspond to P_88. In my opinion, expressions such as “Another line graph similar to Panel A but with different data” don’t communicate this to the reader optimally. Instead, I would propose inserting a sentence about this (similarly as in the alt text of Fig. 4): “Each column represents different conditions with varying pressure values: P12, P50, and P88.” Then, separate descriptions of panels B and C (and similarly E and F and H and I) could be removed.

3.3 SI model reproduced an empirical B. pendula vulnerability curve

So far, we have concentrated on how well the SI model replicated the outcome of the more detailed, physiological embolism spreading model. However, such physiological models are not available for all plant species due to limited availability of data on xylem physiology and anatomy. To demonstrate the potential of the SI model also in these cases, we next fitted the SI model against an empirical B. pendula vulnerability curve. Despite the difference in $P_{50}$ between the physiological model and empirical data, the re-fitted SI model produced a vulnerability curve that closely matched the empirical one (Figure 6). The optimized $\beta$ varied between 0 and 0.11 and increased almost monotonically as a function of $P$ , excluding small noise at the least negative $P$ values and a peak in the beginning of the saturation zone. Both of these were probably artifacts related to the limited number of simulation iterations (100) used for optimizing the $\beta$ value for each $P$ value.

After re-fitting, the SI model allowed approximating the behavior of network properties, for example, prevalence, LCC size, and nonfunctional sap volume, during embolism spreading (Figure S3). This behavior is not directly measurable and cannot be addressed through the empirical vulnerability curve.

Figure 6.

The SI model produces a vulnerability curve that matches empirical observations. A) The empirical vulnerability curve (red line) of B. pendula is from González-Muñoz et al. (Reference González-Muñoz, Sterck, Torres-Ruiz, Petit, Cochard, von and Heinze2018) and obtained as a sigmoidal fit to data. The modeled vulnerability curve (gray) is obtained by optimizing the SI spreading probability ( $\beta$ ) for each $P$ . B) The optimized $\beta$ values used for constructing the SI vulnerability curve.

Two line graphs depict the relationship between the xylem pressure, the PLC extracted from empirical data and simulated with the SI model, and optimized beta.

4. Discussion

4.1 SI model describes embolism spreading based on a single free parameter

In the present article, we have demonstrated that a stochastic model for embolism spreading in plant xylem, the SI model, produces similar results than previous, physiological embolism spreading models. Importantly, in contrast to physiological models that require detailed data on intervessel pit membrane function and anatomy, the SI model depends on a single spreading probability parameter that can be fitted against an empirical vulnerability curve. Therefore, the SI model opens new insights on embolism spreading without requiring detailed information about intervessel pit membrane structure.

The similarity between the physiological and SI models is not limited to the PLC measured at the final state of embolism spreading. Instead, the development of several xylem network properties across time was qualitatively similar in both models. While the earlier saturation in the physiological model and higher variance in the SI model yield rather high normalized RMSE values (see Table 3), the SI model is still able to qualitatively describe the dynamics of the xylem network structure during embolism spreading. While these dynamics can be empirically determined based on time-consuming computed tomography and image analysis (see, e.g., Choat et al. (Reference Choat, Badel, Burlett, Delzon, Cochard and Jansen2016)), they are not accessible based on empirical vulnerability curves. Particularly the shape of the prevalence curve, fast initial growth followed by long saturation, is typical for the SI model and other compartment models without recovery (Pastor-Satorras et al., Reference Pastor-Satorras, Castellano, Van Mieghem and Vespignani2015). Observing this shape of prevalence curve for both the physiological and SI models strongly supports modeling embolism spreading with the SI model.

The most notable difference between the SI and physiological models is the time scale of the saturation of embolism spreading: the physiological model requires around 10-fold less time to saturate than the SI model. However, mapping the simulation time steps to absolute time is not straightforward. Both spreading models consider embolization of a vessel immediate so that all vessels embolized at time $t$ can spread embolism further at time $t + 1$ . However, this is a simplification: in reality, the time required for a vessel to become fully embolized and able to spread the embolism to its neighbors is not constant but varies depending on, for example, the sap pressure, status of the neighboring vessels, and sap lipid concentration. Therefore, we consider the ability of the SI model to reproduce the outcomes of the physiological model qualitatively more important than the difference in the simulation timescales. In particular, we do not recommend using either of the spreading models to estimate the absolute time scale of embolism propagation.

Further, it is important to understand the different meaning of saturation in the two models. In both models, embolism spreading slows down when the number of susceptible vessels decreases, as described above. In the deterministic physiological model, this leads to true saturation: at some point, all vessels that can be embolized at a given $P$ are embolized and the spreading stops. The SI model, on the other hand, would, due to its stochasticity, yield 100% prevalence in an infinitely long simulation. Therefore, at the saturation point, the spreading does not stop but becomes extremely slow.

In our simulations, already relatively low spreading probabilities (around 0.2) yield very high prevalences at the saturation point. Therefore, fitting the model for $P$ values where the full network does not become embolized has yielded very low spreading probability ( $\beta$ ) values (0.01, 0.04, and 0.075 for $P_{12}$ , $P_{50}$ , and $P_{88}$ , respectively) associated with long saturation times. Larger $\beta$ s would guarantee faster saturation but also produce undesirably low final $K_{eff}$ . Another way to shorten the runtime of the SI model would be to record the number of simulation steps the physiological model needed to saturate, $n_{sat}$ , and optimize $\beta$ based on the $K_{eff, SI}(n_{sat})$ yielded by the SI model after $n_{sat}$ steps. However, while this strategy would make the simulations slightly faster, it would weaken the correspondence of the behavior of network properties during the spreading process between the two models. Further, this optimization strategy cannot be used when optimizing the spreading probability against an empirical vulnerability curve as the empirical saturation time is typically unknown.

4.2 Xylem network breakdown contributes to hydraulic failure

In both the physiological and the SI model, the LCC of the xylem network becomes nonfunctional, that is, remains sap-filled but loses its connection to network inlet or outlet and thus stops contributing to sap flow, during embolism spreading. Further, for more negative $P$ s, the final prevalence, or fraction of embolized vessels, is lower than the PLC. In other words, especially at more negative pressures, hydraulic conductivity is lost faster than the xylem sap volume becomes reduced. These results indicate that while our hypothesis about the network breakdown as the main cause of hydraulic failure holds only for the more negative $P$ s, the disruption of the xylem network structure that cuts down the sap pathways is an equally important reason for hydraulic failure as the drying up of the xylem at all pressure values. This observation is supported by earlier results on angiosperms, where leaf water potential and transpiration reduce faster than the fraction of sap-filled vessels in stem xylem (Gauthey et al., Reference Gauthey, Peters, Lòpez, Carins-Murphy, Rodriguez-Dominguez, Tissue and Choat2022), and conifers, where xylem pressures close to $P_{50}$ cause 25% loss in branch tissue water content (Rosner et al., Reference Rosner, Nöbauer and Voggeneder2021). Our observations on the nonfunctional sap volume confirm this result: significant amounts of sap remain in the xylem tissue even when 88% of conductance is lost.

While the final nonfunctional sap volume indicates presence of sap in the xylem also after hydraulic failure in both models, the behavior of nonfunctional sap volume during the spreading process differs between the models. In the physiological model, we observe a nearly monotonic, continuous increase, while in the SI model, the nonfunctional sap volume is initially low and increases only close to the saturation point. This difference reflects the deterministic and stochastic natures of the two models. Let $\mathcal{C}$ be a part of the xylem that is connected to the rest of the network, and thus to inlet and outlet, only through vessel $I$ . Now, embolization of $I$ cuts $\mathcal{C}$ apart from the network, and it becomes a nonfunctional component. In the physiological model, embolism cannot further spread to $\mathcal{C}$ if the BPP between $I$ and its neighbors in $\mathcal{C}$ is more negative than $P$ . Therefore, the total volume of vessels in $\mathcal{C}$ contributes to increase of the nonfunctional sap volume until the end of the simulation. In the SI model, on the other hand, embolism can still spread from $I$ to $\mathcal{C}$ with probability $\beta$ , and the nonfunctional sap volume decreases as vessels in $\mathcal{C}$ become embolized. At the level of the whole network, the simultaneous addition of new nonfunctional components and removal of embolized vessels in existing nonfunctional components keeps the overall nonfunctional sap volume lower than in the physiological model. At very negative $P$ values, BPP values more negative than $P$ become rare and embolism can thus spread to most of the nonfunctional components. Indeed, at more negative $P$ s, the behavior of nonfunctional sap volume in the physiological model approaches that observed in the SI model.

4.3 Network structure and model parameters affect modeled embolism spreading

Historically, intervessel connectivity has been seen as a safety mechanism, following the reasoning of Carlquist (Reference Carlquist1977): in presence of numerous, parallel vessels, an embolism in a single vessel is less probable to cause hydraulic failure in the entire xylem. However, under certain conditions tight intervessel connectivity can rather increase than decrease vulnerability of the xylem (Ewers et al., Reference Ewers, Ewers, Jacobsen and López-Portillo2007; Loepfe et al., Reference Loepfe, Martinez-Vilalta, Piñol and Mencuccini2007; Wheeler et al., Reference Wheeler, Sperry, Hacke and Hoang2005). From the viewpoint of pit membrane structure, this finding is in line with the pit area hypothesis (i.e., large pit membrane areas are more vulnerable to embolism spreading than small areas) but not with the rare pit hypothesis (i.e., embolism resistance depends on a single, exceptionally large pit membrane pore) (Kaack et al., Reference Kaack, Weber, Isasa, Karimi, Li, Pereira and Jansen2021; Wheeler et al., Reference Wheeler, Sperry, Hacke and Hoang2005). In general, this finding indicates a more complex relationship between intervessel networks, sap flow, and embolism spreading in xylem than assumed before. Meanwhile, modern advanced imaging techniques, such as computed tomography (Brodersen et al., Reference Brodersen, Lee, Choat, Jansen, Phillips, Shackel and Matthews2011; Choat et al., Reference Choat, Badel, Burlett, Delzon, Cochard and Jansen2016), have yielded new understanding of the intricacies of xylem architecture. Embolism propagation is a dynamic process that takes place on top of the intervessel network. Thus it is clear that appropriate representation of the complex 3D structure of the xylem will be critical to accurately simulate embolism propagation.

The behavior of prevalence as a function of time depends on network structure and average degree of network nodes (Barthélemy et al., Reference Barthélemy, Barrat, Pastor-Satorras and Vespignani2005; Pastor-Satorras et al., Reference Pastor-Satorras, Castellano, Van Mieghem and Vespignani2015; Vespignani, Reference Vespignani2012). In real-world networks, also network size affects the prevalence. In the present article, we used relatively small networks (100 × 10 × 100 potential vessel elements). We did not investigate the effects of the simulated xylem network size on $P_{50}$ . Therefore, applying the model on networks of a different size should be considered with care.

In the case of physiological models, the final state of embolism spreading is affected both by the underlying xylem network structure and the modeling of intervessel pit membrane function. In our physiological model, we modeled the intervessel pit membrane as a 3D object instead of the more traditional 2D surface model. Importantly, in the 3D model, intervessel pit membrane pores are not 2D openings but channels that comprise multiple pore constrictions. This affects the estimation of BPP across the pit membrane: while most 2D models, including the one used by Mrad et al. (Reference Mrad, Domec, Huang, Lens and Katul2018) and Mrad et al. (Reference Mrad, Johnson, Love and Domec2021), build on the rare pit hypothesis, simulations using the 3D models show that the probability to find exceptionally large pit membrane pores is very low especially in species with thicker pit membranes (Kaack et al., Reference Kaack, Weber, Isasa, Karimi, Li, Pereira and Jansen2021). Therefore, the 2D models may produce artificially high embolism propagation probabilities and overestimate species’ vulnerability to hydraulic failure due to embolism. The intervessel pit membrane of B. pendula is relatively thin (205 nm according to Kaack et al. (Reference Kaack, Weber, Isasa, Karimi, Li, Pereira and Jansen2021)). Therefore, the effect of the pit membrane model selection may not be particularly large in our analysis. However, the difference between the 2D and 3D pit models becomes more prominent when the pit membrane thickness increases. Therefore, it is a notable source of error for many angiosperms (Kaack et al., Reference Kaack, Weber, Isasa, Karimi, Li, Pereira and Jansen2021). Empirical observations on the movement of nanoparticles through pit membranes support the selection of the 3D model (Kaack et al., Reference Kaack, Weber, Isasa, Karimi, Li, Pereira and Jansen2021; Zhang et al., Reference Zhang, Pereira, Kaack, Liu and Jansen2024).

The $P_{50}$ value of our physiological model, −1.28 MPa, was slightly less negative than the empirically measured values for B. pendula, which varied from −1.59 to −1.9 for six vulnerability curves fitted based on the flow-centrifuge measurements (Dulamsuren et al., Reference Dulamsuren, Abilova, Bektayeva, Eldarov, Schuldt, Leuschner and Hauck2019; González-Muñoz et al., Reference González-Muñoz, Sterck, Torres-Ruiz, Petit, Cochard, von and Heinze2018). We consider these values largely comparable because a 10% to 30% error from the true $P_{12}$ , $P_{50}$ , and $P_{88}$ values is acceptable and equal to the typical disagreement in embolism resistance measured by different methods on the same species (Silva, Pereira et al., Reference Silva, Pereira, Kaack, Guan, Pfaff, Trabi and Jansen2024; Silva, Pfaff et al., Reference Silva, Pfaff, Pereira, Miranda and Jansen2024; Yang et al., Reference Yang, Pereira, Peng, Ribeiro, Kaack, Jansen and Tyree2023). However, there are also potential error sources that may explain the difference between the modeled and empirical $P_{50}$ values. First, the embolism resistance measured by González-Muñoz et al. (Reference González-Muñoz, Sterck, Torres-Ruiz, Petit, Cochard, von and Heinze2018) could have been slightly overestimated due to time-based shifts in vulnerability curves generated with a flow-centrifuge (Silva, Pfaff et al., Reference Silva, Pfaff, Pereira, Miranda and Jansen2024) and the temporal buffering effect of water release from embolism on water potential declines and embolism spreading (Hölttä et al., Reference Hölttä, Cochard, Nikinmaa and Mencuccini2009; Vergeynst et al., Reference Vergeynst, Dierick, Bogaerts, Cnudde and Steppe2015). However, the potential buffering effect is minor and needs to be investigated in further detail.

On the other hand, the difference between the modeled and empirical $P_{50}$ values may relate also to the pit membrane model, which yielded for single pit membranes BPP values slightly higher than the empirical $P_{50}$ values of the same species, as suggested also by Kaack et al. (Reference Kaack, Weber, Isasa, Karimi, Li, Pereira and Jansen2021). This indicates that despite its more realistic view of the pit membrane as a 3D object, the model still describes some interactions between liquids, surfactants, and gas nanobubbles at the pit membrane inaccurately. Understanding these interactions is an important direction for future research.

Further, our physiological model selects the seed vessels for embolism spreading in a way that is not entirely realistic. First, a certain amount of embolized vessels is known to exist in intact plants at native pressures (Torres-Ruiz et al., Reference Torres-Ruiz, Cochard, Mencuccini, Delzon and Badel2016), while understanding the origin of these initial embolisms remains a subject for further research (Carmesin et al., Reference Carmesin, Port, Böhringer, Gottschalk, Rasche and Jansen2023). Therefore, selecting a single seed per network component may be too little. Second, while we select the seed vessels at random, in reality these vessels are more common in the most central parts of the xylem (Choat et al., Reference Choat, Badel, Burlett, Delzon, Cochard and Jansen2016). However, underestimating the number of seed vessels cannot explain the observed difference between the empirical and observed values, since it should make the modeled $P_{50}$ value more, not less, negative. Further, the effect of ignoring the distribution of seed vessel locations is probably negligible, as in our model the properties of seed vessels, most importantly their degree (i.e., number of neighbors) that affects the SI spreading process, are independent of the location in the xylem.

In addition to the difference between the $P_{50}$ values, the output of our physiological model differed from empirical data in terms of the vulnerability curve shape. Our model produced $P_{12}$ and $P_{88}$ values closer to each other than observed empirically (−1.26 MPa and −1.29 MPa vs. −1.62 MPa and −2.00 MPa in González-Muñoz et al. (Reference González-Muñoz, Sterck, Torres-Ruiz, Petit, Cochard, von and Heinze2018)). Thus, in our simulations PLC increased as a function of $P$ faster than what has been observed empirically, and the transition from normal xylem function to hydraulic failure happened over an unrealistically narrow range of pressures. Since our investigated pressure range was particularly dense around the transition point (for details, see Section 2.3.3), the steep increase in PLC is a feature of the physiological model itself. Most probably, this discrepancy is related to the fact that the model considers embolism events immediate so that a vessel can spread embolism further right after becoming embolized. In reality, embolism spreads from a newly embolized vessel only after a certain gas concentration is reached (Pereira et al., Reference Pereira, Kaack, Guan, Silva, Miranda, Pires, Ribeiro, Schenk and Jansen2023; Silva et al., Reference Silva, Bujnowski, Pereira, Miranda, Schenk and Jansen2025), and ignoring this delay may make the vulnerability curve more step-like. In future, this problem can be addressed by adding to the model a stochastic waiting time, during which a newly embolized vessel cannot yet spread embolism, and fine-tuning the related compartment model (for details, see Section 4.4).

Finally, our model concentrated on embolism spreading in a network formed by a single type of conduits, vessels. In B. pendula, imperforate tracheary elements (ITEs) are unlikely to be conductive (as Sano et al. (Reference Sano, Morris, Shimada, De Craene LP and Jansen2011) observed for a related species, B. japonica). However, they may still affect the resistance against embolism propagation through their mechanistic support role, and excluding them from the model is thus a potential, although minor, source of error. Further, including ITEs in the model would yield a more realistic picture of embolism spreading for the many angiosperm species with ITEs as a second type of conduits in addition to vessels (Esteban et al., Reference Esteban, de, Gasson, Garcıa-Iruela, Garcıa-Fernández and Garcıa-Esteban2024; Olson, Reference Olson2023). To this end, a promising approach would be a multilayer network model (Kivelä et al., Reference Kivelä, Arenas, Barthelemy, Gleeson, Moreno and Porter2014), where networks of different conducting cells, including vessels and ITEs, would form layers of their own while the interconnectivity between conduits of different types would be represented by the interlayer links.

4.4 Which compartment model is best suited for modeling embolism spreading?

The SI model belongs to the family of compartment models. These models divide individuals (here, xylem vessels) into several groups, or compartments, and monitor the transitions between these groups (Hethcote, Reference Hethcote2000; Pastor-Satorras et al., Reference Pastor-Satorras, Castellano, Van Mieghem and Vespignani2015). As the simplest member of the family, the SI model comprises only two compartments: the susceptible (S) (i.e., non-embolized) and infected (I) (i.e., embolized) individuals. More complicated compartment models include, for example, recovery from infection by returning to the susceptible compartment (the SIS model) or removal of infected individuals through immunity or death (the SIR model) (Pastor-Satorras et al., Reference Pastor-Satorras, Castellano, Van Mieghem and Vespignani2015).

The SI model is a strong simplification in the case of many infectious diseases (Pastor-Satorras et al., Reference Pastor-Satorras, Castellano, Van Mieghem and Vespignani2015). However, it describes embolism spreading more naturalistically than other compartment models. For example, the SIS model would suggest recovery from embolism through refilling of embolized vessels. However, such refilling is rarely observed in woody plants under negative pressure (Choat et al., Reference Choat, Brodersen and McElrone2015; Choat et al., Reference Choat, Nolf, Lopez, Peters, Carins-Murphy, Creek and Brodribb2019; Gauthey et al., Reference Gauthey, Peters, Lòpez, Carins-Murphy, Rodriguez-Dominguez, Tissue and Choat2022; Knipfer et al., Reference Knipfer, Brodersen, Zedan, Kluepfel, McElrone and Phillips2015), and even in non-woody plants these mechanisms take place at time scales longer than those considered in our spreading simulations (Stewart et al., Reference Stewart, Allen, Polutchko, Ocheltree and Gleason2025). In B. pendula, refilling takes place under seasonal positive pressures, but it remains unclear if this is possible under normal, negative xylem pressures (Choat et al., Reference Choat, Brodribb, Brodersen, Duursma, López and Medlyn2018).

On the other hand, as discussed in Section 4.3, a newly embolized vessel can spread embolism further only after a waiting time. Therefore, the model could be further improved by adding an exposed (E) compartment for vessels that are already embolized but do not yet spread the embolism. In this hypothetical SEI model, the E–I transition could happen at a constant rate based on empirical observations. On the other hand, if the S-E transition was defined as a snap-off event yielding to an increased number of gas nanobubbles in the exposed vessel, the E–I transition could be also stochastic and defined by the probability of nanobubble expansion (Ingram et al., Reference Ingram, Reischl, Vesala and Vehkamäki2024). Another potential extension of the SI model would be modeling several, possibly interacting, spreading processes on the same xylem network. Such model could describe, for example, the simultaneous spreading of embolism and fungi in the xylem.

Earlier, Roth-Nebelsick (Reference Roth-Nebelsick2019) has modeled embolism spreading with the SIR model, where the removed (R) compartment consisted of areas that can no more spread embolism due to lack of susceptible neighbors. Torres-Ruiz et al. (Reference Torres-Ruiz, Cochard, Mencuccini, Delzon and Badel2016), on the other hand, has modeled embolism spreading between tracheids using an SI model variant that allowed both spontaneous embolism and spreading between orthogonally adjacent tracheids. While these models were able to reproduce empirical results, they considered the xylem as a continuous medium, ignoring the intervessel connectivity. However, the connectivity structure plays an important role in xylem hydraulics (Avila et al., Reference Avila, Kane, Batz, Trabi, Damatta, Jansen and McAdam2023; Cai et al., Reference Cai, Li, Zhang, Zhang and Tyree2014; Guan et al., Reference Guan, Pereira, McAdam, Cao and Jansen2021; Loepfe et al., Reference Loepfe, Martinez-Vilalta, Piñol and Mencuccini2007) and in spreading processes in general (Johnson, Reference Johnson2024; Newman, Reference Newman2002; Pastor-Satorras et al., Reference Pastor-Satorras, Castellano, Van Mieghem and Vespignani2015). Our approach, running the SI model on top of simulated xylem networks, combines the small number of free parameters typical for compartment models with the possibility to investigate how intervessel connectivity affects embolism spreading and to observe the changes in the xylem network structure during the spreading process.

4.5 Embolism spreading is a directed percolation phenomenon

In addition to exploring the effects of intervessel connectivity on xylem hydraulics, running spreading models on top of the intervessel network allows investigating embolism spreading with tools adopted from network science and applied to investigate a diverse set of real-world systems (Newman, Reference Newman2003). From this viewpoint, embolism spreading is an instance of a well-known concept of theoretical physics, directed (bond) percolation, that is widely applied to study inherently directed phenomena such as filtering, or percolation, of fluids through a porous media (Badie-Modiri et al., Reference Badie-Modiri, Rizi, Karsai and Kivelä2022; Hinrichsen, Reference Hinrichsen2000). In network science, a theoretical mapping has been established between directed percolation and SI spreading processes on top of networks (Badie-Modiri et al., Reference Badie-Modiri, Rizi, Karsai and Kivelä2022; Kenah and Robins, Reference Kenah and Robins2007; Newman, Reference Newman2002; Pastor-Satorras et al., Reference Pastor-Satorras, Castellano, Van Mieghem and Vespignani2015).

Directed percolation processes are characterized by a sudden transition between permeable and impermeable phases (Badie-Modiri et al., Reference Badie-Modiri, Rizi, Karsai and Kivelä2022; Hinrichsen, Reference Hinrichsen2000). These transitions are observed in terms of two parameters: as the value of the free control parameter changes, the system transitions between regions where an observable order parameter gets either zero or non-zero values (Badie-Modiri et al., Reference Badie-Modiri, Rizi, Karsai and Kivelä2022; Hinrichsen, Reference Hinrichsen2000). In network science, the two phases correspond to connected and fragmented network structures, and the percolation transition is observed when a network breaks down due to link removal (Newman, Reference Newman2003). Here, the control parameter is the fraction of removed links, while the most commonly observed order parameter is the size of the giant component, defined as a component that spans a finite fraction of the network when the network size grows to infinity. In empirical studies, the size of the LCC approximates the giant component size. As predicted by the percolation theory, most networks are relatively robust against link removal until a critical fraction of links is achieved. After this, the network quickly breaks apart, as indicated by a sudden drop in the LCC size. The exact value of the critical fraction of removed links, known as the percolation threshold, depends on the degree distribution of the network (Newman, Reference Newman2003).

When considering the SI spreading process as directed percolation, the control and order parameters correspond to the spreading probability and fraction of infected nodes, correspondingly. The percolation is directed because of the causal nature of infection spreading: while the infection can spread along network links in all spatial directions, the spreading along an additional time axis is limited, since nodes can infect their neighbors only after they themselves have been infected. Because of the lack of recovery mechanisms in the SI model, any non-zero spreading probability will eventually lead to a fully infected network. Therefore, the two states, between which the transition takes place, are not defined by the final prevalence but the prevalence at the saturation phase, where spreading slows down significantly (see Figures 4 and 5). If this saturation happens only when prevalence is approximately 100%, the entire network will become infected in a finite time also when the network size grows to infinity, leading to non-zero order parameter values. On the other hand, if the saturation happens at a lower prevalence value, the time required to infect the entire network grows to infinity together with network size.

In embolism spreading, we see the sigmoidal shape characteristic for the percolation transition in the vulnerability curves that visualize the PLC as a function of $P$ (see Figure 3A and, e.g., Avila et al. (Reference Avila, Guan, Kane, Cardoso, Batz, DaMatta and McAdam2022); González-Muñoz et al. (Reference González-Muñoz, Sterck, Torres-Ruiz, Petit, Cochard, von and Heinze2018)). This is not surprising, as $K_{eff}$ depends on the percolation order parameter, the LCC size. This shape indicates that the plant xylem is relatively resistant to embolism until a certain threshold is reached, after which the embolism quickly spans the whole vessel network, leading to a percolation-like collapse of network structure (Avila et al., Reference Avila, Guan, Kane, Cardoso, Batz, DaMatta and McAdam2022; Delzon and Cochard, Reference Delzon and Cochard2014). There are two explanations for the sigmoidal shape of the vulnerability curve. First, in all spreading processes, the spreading speed increases with the number of infected individuals since, in the early stages, the number of susceptible individuals (here, sap-filled vessels) with infected neighbors increases faster than the number of infected individuals (here, embolized vessels) itself (Pastor-Satorras et al., Reference Pastor-Satorras, Castellano, Van Mieghem and Vespignani2015). Second, in plants specifically, above the embolism formation threshold, embolism propagation depends less strictly on the pressure difference across intervessel pit membranes, possibly because of the increased role of gas diffusion in embolism spreading at more negative $P$ values (Avila et al., Reference Avila, Guan, Kane, Cardoso, Batz, DaMatta and McAdam2022; Silva, Pereira et al., Reference Silva, Pereira, Kaack, Guan, Pfaff, Trabi and Jansen2024). This decreases the effect of pit safety mechanisms and thus makes the xylem more vulnerable to embolism spreading (Avila et al., Reference Avila, Guan, Kane, Cardoso, Batz, DaMatta and McAdam2022).

The role of gas diffusion in speeding up embolism spreading is two-fold. On one hand, after becoming embolized, a vessel needs to reach the atmospheric gas concentration before spreading embolism to neighboring vessels (Pereira et al., Reference Pereira, Kaack, Guan, Silva, Miranda, Pires, Ribeiro, Schenk and Jansen2023; Silva et al., Reference Silva, Bujnowski, Pereira, Miranda, Schenk and Jansen2025). This typically happens through diffusion both radially and from neighboring sap-filled vessels, which, as discussed in Section 4.3, may act as a temporary buffer mechanism against embolism in the neighboring vessels (Silva et al., Reference Silva, Bujnowski, Pereira, Miranda, Schenk and Jansen2025). However, when the fraction of embolized vessels exceeds the percolation threshold, the high availability of air in the xylem speeds up diffusion to newly embolized vessels, thus weakening this buffer mechanism (Avila et al., Reference Avila, Kane, Batz, Trabi, Damatta, Jansen and McAdam2023; Silva et al., Reference Silva, Bujnowski, Pereira, Miranda, Schenk and Jansen2025; Silva, Pereira et al., Reference Silva, Pereira, Kaack, Guan, Pfaff, Trabi and Jansen2024). On the other hand, after the embolized vessels have reached the atmospheric gas concentration, they become air sources themselves and can spread embolism both through snap-off events and diffusion (Silva, Pfaff et al., Reference Silva, Pfaff, Pereira, Miranda and Jansen2024). Since the detailed mechanisms related to diffusion in embolism spreading are still not fully known, they are not included in our physiological model. However, fitting the SI model against empirical data can account for the increased role of diffusion through an increased spreading probability, $\beta$ , above the embolism formation threshold. Interestingly, central evidence for the role of diffusion comes from the increase of PLC as a function of time when $P$ is kept constant. While the SI model does not explicitly model diffusion, similar increase of PLC over time is seen due to the stochastic nature of the model.

Defining the exact mapping between embolism spreading and directed percolation remains as a subject for future research. However, our results on the ability of the SI model to describe embolism spreading in xylem, together with the known theoretical equivalence between SI spreading and directed percolation, demonstrate that embolism spreading belongs to the universality class of directed percolation. Mrad et al. (Reference Mrad, Johnson, Love and Domec2021) has already proposed applying percolation theory to understand the connection between the average connectivity between vessels and embolism threshold. Finding the exact mapping to directed percolation will open further application possibilities for percolation theory, including the analytical and numerical results on the size distribution of embolized components and behavior of prevalence across time.

5. Conclusion

We have demonstrated that a stochastic model for embolism spreading, the SI model, produces vulnerability curves comparable to those produced by a physiological model or extracted from empirical data. The succesful fitting of the SI model against the empirical vulnerability curve of B. pendula shows that the model can be used to investigate embolism spreading dynamics also in species, for which detailed data on xylem physiology are not available. Running spreading models on top of a simulated xylem network allows exploring the effects of intervessel connectivity on embolism spreading and monitoring xylem network properties during the spreading process. In particular, we observed that during the spreading, the largest connected component of the intervessel network becomes nonfunctional and does not further contribute to sap flow. This indicates that conductance loss is caused by both direct embolism and groups of non-embolized vessels being cut off from the functional xylem network. Combining the SI spreading model with intervessel connectivity structure also sets embolism spreading into a broader context of network phenomena as an instance of directed percolation.

Acknowledgements

We thank Jari Saramäki and Mikko Kivelä for fruitful discussions on the SI spreading model. Further, we gratefully acknowledge the computational resources provided by the Aalto Science-IT project.

Supplementary material

The supplementary material for this article can be found at http://doi.org/10.1017/nws.2026.10035.

Competing interests

None.

Funding statement

This work was supported by the Research Council Finland (O.K., P.K., A.L., and AM.L., decision number 342934), (M.H. and A.L., decision numbers 355142 and 359340).

Data availability statement

All analysis of the present study are based on simulated data. Python code for simulating the data and performing analyses is available at GitHub: https://github.com/onerva-korhonen/tree-percolation

Ethical standards

The research meets all ethical guidelines, including adherence to the legal requirements of the study country.

Author contributions

Conceptualization: O.K.; S.J.; AM.L. Data curation: O.K. Formal analysis: O.K. Funding acquisition: A.L., AM.L. Investigation: O.K. Methodology: O.K., S.J., L.M.S., P.K. Project administration: O.K., AM.L. Resources: S.J., A.L., M.H. Software: O.K., P.K. Supervision: AM.L. Validation: O.K. Visualization: O.K. Writing - original draft: O.K. Writing - review & editing: all authors. All authors approved the final submitted draft.

References

Adams, H. D., Zeppel, M. J., Anderegg, W. R., Hartmann, H., Landhäusser, S. M., Tissue, D. T. … Allen, C. D. (2017). A multi-species synthesis of physiological mechanisms in drought-induced tree mortality. Nature Ecology & Evolution, 1(9), 1285–1291.10.1038/s41559-017-0248-xCrossRef Google Scholar PubMed

Alber, M., Petit, G., & Sellin, A. (2019). Does elevated air humidity modify hydraulically relevant anatomical traits of wood in Betula pendula? Trees, 33, 1361–1371.10.1007/s00468-019-01863-0CrossRef Google Scholar

Anderegg, W. R., Klein, T., Bartlett, M., Sack, L., Pellegrini, A. F., Choat, B. … Jansen, S. (2016). Meta-analysis reveals that hydraulic traits explain cross-species patterns of drought-induced tree mortality across the globe. Proceedings of the National Academy of Sciences, 113(18), 5024–5029.10.1073/pnas.1525678113CrossRef Google Scholar PubMed

Avila, R. T., Guan, X., Kane, C. N., Cardoso, A. A., Batz, T. A., DaMatta, F. M. … McAdam, S. A. (2022). Xylem embolism spread is largely prevented by interconduit pit membranes until the majority of conduits are gas-filled. Plant, Cell & Environment, 45(4), 1204–1215.10.1111/pce.14253CrossRef Google Scholar PubMed

Avila, R. T., Kane, C. N., Batz, T. A., Trabi, C., Damatta, F. M., Jansen, S. … McAdam, S. A. (2023). The relative area of vessels in xylem correlates with stem embolism resistance within and between genera. Tree Physiology, 43(1), 75–87.10.1093/treephys/tpac110CrossRef Google Scholar PubMed

Badie-Modiri, A., Rizi, A. K., Karsai, M., & Kivelä, M. (2022). Directed percolation in random temporal network models with heterogeneities. Physical Review E, 105(5), 054313.10.1103/PhysRevE.105.054313CrossRef Google Scholar PubMed

Barrat, A., Cattuto, C., Kivelä, M., Lehmann, S., & Saramäki, J. (2021). Effect of manual and digital contact tracing on COVID-19 outbreaks: a study on empirical contact data. Journal of the Royal Society Interface, 18(178), 20201000.10.1098/rsif.2020.1000CrossRef Google Scholar

Barthélemy, M., Barrat, A., Pastor-Satorras, R., & Vespignani, A. (2005). Dynamical patterns of epidemic outbreaks in complex heterogeneous networks. Journal of Theoretical Biology, 235(2), 275–288.10.1016/j.jtbi.2005.01.011CrossRef Google Scholar PubMed

Bhat, K., & Kärkkäinen, M. (1981). Variation in structure and selected properties of Finnish birch wood. IV. Silva Fennica, 15(1), 10–17.Google Scholar

Bosshard, H. H., & Kucera, L. (1973). Die dreidimensionale Strukturanalyse des holzes—Erste Mitteilung: Die Vernetzung des Gefäßsystems in Fagus sylvatica L. European Journal of Wood and Wood Products, 31(11), 437–445.10.1007/BF02613953CrossRef Google Scholar

Brodersen, C. R., Lee, E. F., Choat, B., Jansen, S., Phillips, R. J., Shackel, K. A. … Matthews, M. A. (2011). Automated analysis of three-dimensional xylem networks using high-resolution computed tomography. New Phytologist, 191(4), 1168–1179.10.1111/j.1469-8137.2011.03754.xCrossRef Google Scholar PubMed

Burggraaf, P. (1972). Some observations on the course of the vessels in the wood of Fraxinus excelsior L. Acta Botanica Neerlandica, 21(1), 32–47.10.1111/j.1438-8677.1972.tb00744.xCrossRef Google Scholar

Cai, J., Li, S., Zhang, H., Zhang, S., & Tyree, M. T. (2014). Recalcitrant vulnerability curves: methods of analysis and the concept of fibre bridges for enhanced cavitation resistance. Plant, Cell & Environment, 37(1), 35–44.10.1111/pce.12120CrossRef Google Scholar PubMed

Carlquist, S. (1977). Ecological factors in wood evolution: A floristic approach. American Journal of Botany, 64(7), 887–896.10.1002/j.1537-2197.1977.tb11932.xCrossRef Google Scholar

Carmesin, C. F., Port, F., Böhringer, S., Gottschalk, K. E., Rasche, V., & Jansen, S. (2023). Ageing-induced shrinkage of intervessel pit membranes in xylem of Clematis vitalba modifies its mechanical properties as revealed by atomic force microscopy. Frontiers in Plant Science, 14, 1002711.10.3389/fpls.2023.1002711CrossRef Google Scholar PubMed

Choat, B., Badel, E., Burlett, R., Delzon, S., Cochard, H., & Jansen, S. (2016). Noninvasive measurement of vulnerability to drought-induced embolism by X-ray microtomography. Plant Physiology, 170(1), 273–282.10.1104/pp.15.00732CrossRef Google Scholar PubMed

Choat, B., Brodersen, C. R., & McElrone, A. J. (2015). Synchrotron X-ray microtomography of xylem embolism in Sequoia sempervirens saplings during cycles of drought and recovery. New Phytologist, 205(3), 1095–1105.10.1111/nph.13110CrossRef Google Scholar PubMed

Choat, B., Brodribb, T. J., Brodersen, C. R., Duursma, R. A., López, R., & Medlyn, B. E. (2018). Triggers of tree mortality under drought. Nature, 558(7711), 531–539.10.1038/s41586-018-0240-xCrossRef Google Scholar PubMed

Choat, B., Nolf, M., Lopez, R., Peters, J. M., Carins-Murphy, M. R., Creek, D. … Brodribb, T. J. (2019). Non-invasive imaging shows no evidence of embolism repair after drought in tree species of two genera. Tree Physiology, 39(1), 113–121.10.1093/treephys/tpy093CrossRef Google Scholar PubMed

Dagan, Z., Weinbaum, S., & Pfeffer, R. (1982). An infinite-series solution for the creeping motion through an orifice of finite length. Journal of Fluid Mechanics, 115, 505–523.10.1017/S0022112082000883CrossRef Google Scholar

Delzon, S., & Cochard, H. (2014). Recent advances in tree hydraulics highlight the ecological significance of the hydraulic safety margin. New Phytologist, 203(2), 355–358.10.1111/nph.12798CrossRef Google Scholar PubMed

Dixon, H. H., & Joly, J. (1895). XII. On the ascent of sap. Philosophical Transactions of the Royal Society of London. (B.), 186, 563–576.10.1098/rstb.1895.0012CrossRef Google Scholar

Duanmu, J. L., & Chai, W. K. (2025). Modelling innovation adoption spreading in complex networks. Applied Network Science, 10(1), 10.10.1007/s41109-025-00698-8CrossRef Google Scholar

Dulamsuren, C., Abilova, S. B., Bektayeva, M., Eldarov, M., Schuldt, B., Leuschner, C. … Hauck, M. (2019). Hydraulic architecture and vulnerability to drought-induced embolism in southern boreal tree species of Inner Asia. Tree Physiology, 39(3), 463–473.10.1093/treephys/tpy116CrossRef Google Scholar PubMed

Emory, S. F. (1989). Principles of integrity-testing hydrophilic microporous membrane filters, part I. Pharmaceutical Technology, 13(9), 68–77.Google Scholar

Esteban, L. G., de, P. P., Gasson, P., Garcıa-Iruela, A., Garcıa-Fernández, F., & Garcıa-Esteban, L. (2024). Hardwoods: Anatomy and functionality of their elements—a short review. Forests, 15(7), 1162.10.3390/f15071162CrossRef Google Scholar

Ewers, F. W., Ewers, J. M., Jacobsen, A. L., & López-Portillo, J. (2007). Vessel redundancy: Modeling safety in numbers. IAWA Journal, 28(4), 373–388.10.1163/22941932-90001650CrossRef Google Scholar

Gauthey, A., Peters, J. M., Lòpez, R., Carins-Murphy, M. R., Rodriguez-Dominguez, C. M., Tissue, D. T. … Choat, B. (2022). Mechanisms of xylem hydraulic recovery after drought in Eucalyptus saligna . Plant, Cell & Environment, 45(4), 1216–1228.10.1111/pce.14265CrossRef Google Scholar PubMed

González-Muñoz, N., Sterck, F., Torres-Ruiz, J. M., Petit, G., Cochard, H., von, A. G. … Heinze, B. (2018). Quantifying in situ phenotypic variability in the hydraulic properties of four tree species across their distribution range in Europe. PLoS One, 13(5), e0196075.10.1371/journal.pone.0196075CrossRef Google Scholar PubMed

Gostick, J., Aghighi, M., Hinebaugh, J., Tranter, T., Hoeh, M. A., Day, H. … Putz, A. (2016). OpenPNM: A pore network modeling package. Computing in Science & Engineering, 18(4), 60–74.10.1109/MCSE.2016.49CrossRef Google Scholar

Guan, X., Pereira, L., McAdam, S. A., Cao, K. F., & Jansen, S. (2021). No gas source, no problem: Proximity to pre-existing embolism and segmentation affect embolism spreading in angiosperm xylem by gas diffusion. Plant, Cell & Environment, 44(5), 1329–1345.10.1111/pce.14016CrossRef Google Scholar PubMed

Held, M., Jyske, T., & Lintunen, A. (2025). Coordinated_scaling_of_conduits_and_pits (Gitlab 4198) [Data set]. Gitlab. https://version.helsinki.fi/magdahel/phd-publication-ii Google Scholar

Held, M., Jyske, T., Lintunen, A., & Sevanto, S. (2026). Conduit and pit dimensions scale in a coordinated way from the treetop to coarse roots in three boreal tree species. Tree Physiology, 46(2), tpag003.10.1093/treephys/tpag003CrossRef Google Scholar

Hethcote, H. W. (2000). The mathematics of infectious diseases. SIAM Review, 42(4), 599–653.10.1137/S0036144500371907CrossRef Google Scholar

Hinrichsen, H. (2000). Non-equilibrium critical phenomena and phase transitions into absorbing states. Advances in Physics, 49(7), 815–958.10.1080/00018730050198152CrossRef Google Scholar

Hölttä, T., Cochard, H., Nikinmaa, E., & Mencuccini, M. (2009). Capacitive effect of cavitation in xylem conduits: results from a dynamic model. Plant, Cell & Environment, 32(1), 10–21.10.1111/j.1365-3040.2008.01894.xCrossRef Google Scholar PubMed

Ingram, S., Reischl, B., Vesala, T., & Vehkamäki, H. (2024). Ruptures of mixed lipid monolayers under tension and supercooling: implications for nanobubbles in plants. Nanoscale Advances, 6(15), 3775–3784.10.1039/D4NA00316KCrossRef Google Scholar PubMed

Johnson, S. (2024). Epidemic modelling requires knowledge of the social network. Journal of Physics: Complexity, 5(1), 01LT01.Google Scholar

Kaack, L., Altaner, C. M., Carmesin, C., Diaz, A., Holler, M., Kranz, C. … Donaldson, L. A. (2019). Function and three-dimensional structure of intervessel pit membranes in angiosperms: A review. IAWA Journal, 40(4), 673–702.10.1163/22941932-40190259CrossRef Google Scholar

Kaack, L., Weber, M., Isasa, E., Karimi, Z., Li, S., Pereira, L. … Jansen, S. (2021). Pore constrictions in intervessel pit membranes provide a mechanistic explanation for xylem embolism resistance in angiosperms. New Phytologist, 230(5), 1829–1843.10.1111/nph.17282CrossRef Google Scholar PubMed

Karimi, Z. (2014). Systematic, Evolutionary and Functional Anatomy of Wood and Leaves in Betulaceae . PhD dissertation, Germany: Universität Ulm, Fakultät für Naturwissenschaften, Institut für Systematische Botanik und Ökologie.Google Scholar

Kenah, E., & Robins, J. M. (2007). Second look at the spread of epidemics on networks. Physical Review E—Statistical, Nonlinear, and Soft Matter Physics, 76(3), 036113.10.1103/PhysRevE.76.036113CrossRef Google Scholar PubMed

Kivelä, M., Arenas, A., Barthelemy, M., Gleeson, J. P., Moreno, Y., & Porter, M. A. (2014). Multilayer networks. Journal of Complex Networks, 2(3), 203–271.10.1093/comnet/cnu016CrossRef Google Scholar

Knipfer, T., Brodersen, C. R., Zedan, A., Kluepfel, D. A., McElrone, A. J., & Phillips, N. (2015). Patterns of drought-induced embolism formation and spread in living walnut saplings visualized using X-ray microtomography. Tree Physiology, 35(7), 744–755.10.1093/treephys/tpv040CrossRef Google Scholar PubMed

Lintunen, A., & Kalliokoski, T. (2010). The effect of tree architecture on conduit diameter and frequency from small distal roots to branch tips in Betula pendula, Picea abies and Pinus sylvestris . Tree Physiology, 30(11), 1433–1447.10.1093/treephys/tpq085CrossRef Google Scholar PubMed

Loepfe, L., Martinez-Vilalta, J., Piñol, J., & Mencuccini, M. (2007). The relevance of xylem network structure for plant hydraulic efficiency and safety. Journal of Theoretical Biology, 247(4), 788–803.10.1016/j.jtbi.2007.03.036CrossRef Google Scholar PubMed

Lu, D., Smith-Martin, C. M., Muscarella, R., Uriarte, M., & Zheng, T. (2025). A spatio-temporal model of embolism propagation in leaf vein networks. AoB PLANTS, 17(4), plaf020.10.1093/aobpla/plaf020CrossRef Google Scholar

Mantova, M., Herbette, S., Cochard, H., & Torres-Ruiz, J. M. (2022). Hydraulic failure and tree mortality: From correlation to causation. Trends in Plant Science, 27(4), 335–345.10.1016/j.tplants.2021.10.003CrossRef Google Scholar PubMed

McDowell, N., Pockman, W. T., Allen, C. D., Breshears, D. D., Cobb, N., Kolb, T. … Williams, D. G. (2008). Mechanisms of plant survival and mortality during drought: Why do some plants survive while others succumb to drought? New Phytologist, 178(4), 719–739.10.1111/j.1469-8137.2008.02436.xCrossRef Google Scholar PubMed

Mrad, A., Domec, J. C., Huang, C. W., Lens, F., & Katul, G. (2018). A network model links wood anatomy to xylem tissue hydraulic behaviour and vulnerability to cavitation. Plant, Cell & Environment, 41(12), 2718–2730.10.1111/pce.13415CrossRef Google Scholar PubMed

Mrad, A., Johnson, D. M., Love, D. M., & Domec, J. C. (2021). The roles of conduit redundancy and connectivity in xylem hydraulic functions. New Phytologist, 231(3), 996–1007.10.1111/nph.17429CrossRef Google Scholar PubMed

Nekovee, M., Moreno, Y., Bianconi, G., & Marsili, M. (2007). Theory of rumour spreading in complex social networks. Physica A: Statistical Mechanics and its Applications, 374(1), 457–470.10.1016/j.physa.2006.07.017CrossRef Google Scholar

Newman, M. E. (2002). Spread of epidemic disease on networks. Physical Review E, 66(1), 016128.10.1103/PhysRevE.66.016128CrossRef Google Scholar PubMed

Newman, M. E. (2003). The structure and function of complex networks. SIAM Review, 45(2), 167–256.10.1137/S003614450342480CrossRef Google Scholar

Olson, M. E. (2023). Imperforate tracheary element classification for studies of xylem structure–function relations. IAWA Journal, 44(3-4), 439–464.10.1163/22941932-bja10125CrossRef Google Scholar

Pastor-Satorras, R., Castellano, C., Van Mieghem, P., & Vespignani, A. (2015). Epidemic processes in complex networks. Reviews of Modern Physics, 87(3), 925–979.10.1103/RevModPhys.87.925CrossRef Google Scholar

Pereira, L., Kaack, L., Guan, X., Silva, L. M., Miranda, M. T., Pires, G. S., Ribeiro, R. V., Schenk, H. J., & Jansen, S. (2023). Angiosperms follow a convex trade-off to optimize hydraulic safety and efficiency. New Phytologist, 240(5), 1788–1801.10.1111/nph.19253CrossRef Google Scholar PubMed

Raponi, S., Khalifa, Z., Oligeri, G., & Di Pietro, R. (2022). Fake news propagation: A review of epidemic models, datasets, and insights. ACM Transactions on the Web (TWEB), 16(3), 1–34.10.1145/3522756CrossRef Google Scholar

Rizi, A. K., Faqeeh, A., Badie-Modiri, A., & Kivelä, M. (2022). Epidemic spreading and digital contact tracing: Effects of heterogeneous mixing and quarantine failures. Physical Review E, 105(4), 044313.10.1103/PhysRevE.105.044313CrossRef Google Scholar PubMed

Rosner, S., Nöbauer, S., & Voggeneder, K. (2021). Ready for screening: fast assessable hydraulic and anatomical proxies for vulnerability to cavitation of young conifer sapwood. Forests, 12(8), 1104.10.3390/f12081104CrossRef Google Scholar

Roth-Nebelsick, A. (2019). It’s contagious: Calculation and analysis of xylem vulnerability to embolism by a mechanistic approach based on epidemic modeling. Trees, 33(5), 1519–1533.10.1007/s00468-019-01891-wCrossRef Google Scholar

Sano, Y., Morris, H., Shimada, H., De Craene LP, R., & Jansen, S. (2011). Anatomical features associated with water transport in imperforate tracheary elements of vessel-bearing angiosperms. Annals of Botany, 107(6), 953–964.10.1093/aob/mcr042CrossRef Google Scholar PubMed

Schenk, H. J., Steppe, K., & Jansen, S. (2015). Nanobubbles: A new paradigm for air-seeding in xylem. Trends in Plant Science, 20(4), 199–205.10.1016/j.tplants.2015.01.008CrossRef Google Scholar PubMed

Sevanto, S., Mcdowell, N. G., Dickman, L. T., Pangle, R., & Pockman, W. T. (2014). How do trees die? A test of the hydraulic failure and carbon starvation hypotheses. Plant, Cell & Environment, 37(1), 153–161.10.1111/pce.12141CrossRef Google Scholar

Shah, D., & Zaman, T. (2011). Rumors in a network: Who’s the culprit? IEEE Transactions on Information Theory, 57(8), 5163–5181.10.1109/TIT.2011.2158885CrossRef Google Scholar

Silva, L. M, Bujnowski, B., Pereira, L., Miranda, M., Schenk, H., & Jansen, S. (2025). Gas diffusion kinetics in relation to embolism formation and propagation in angiosperm xylem: A mini-review. Acta Horticulturae, 1419, 123–134.10.17660/ActaHortic.2025.1419.16CrossRef Google Scholar

Silva, L. M., Pereira, L., Kaack, L., Guan, X., Pfaff, J., Trabi, C. L., & Jansen, S. (2024). The potential link between gas diffusion and embolism spread in angiosperm xylem: Evidence from flow‐centrifuge experiments and modelling. Plant, Cell & Environment, 47(12), 4977–4991.10.1111/pce.15084CrossRef Google Scholar PubMed

Silva, L. M., Pfaff, J., Pereira, L., Miranda, M. T., & Jansen, S. (2024). Embolism propagation does not rely on pressure only: Time-based shifts in xylem vulnerability curves of angiosperms determine the accuracy of the flow-centrifuge method. Tree Physiology, tpae131.Google Scholar

Sperry, J. S., & Hacke, U. G. (2004). Analysis of circular bordered pit function I. Angiosperm vessels with homogenous pit membranes. American Journal of Botany, 91(3), 369–385.10.3732/ajb.91.3.369CrossRef Google Scholar PubMed

Stewart, J. J., Allen, B. S., Polutchko, S. K., Ocheltree, T. W., & Gleason, S. M. (2025). Xylem embolism refilling revealed in stems of a weedy grass. Proceedings of the National Academy of Sciences, 122(13), e2420618122.10.1073/pnas.2420618122CrossRef Google Scholar PubMed

Strati, S., Patiño, S., Slidders, C., Cundall, E. P., & Mencuccini, M. (2003). Development and recovery from winter embolism in silver birch: Seasonal patterns and relationships with the phenological cycle in oceanic Scotland. Tree Physiology, 23(10), 663–673.10.1093/treephys/23.10.663CrossRef Google Scholar PubMed

Strickland, C., Dangelmayr, G., Shipman, P. D., Kumar, S., & Stohlgren, T. J. (2015). Network spread of invasive species and infectious diseases. Ecological Modelling, 309, 1–9.10.1016/j.ecolmodel.2015.04.010CrossRef Google Scholar

Tavares, J. V., Oliveira, R. S., Mencuccini, M., Signori-Müller, C., Pereira, L., Diniz, F. C. … Galbraith, D. R. (2023). Basin-wide variation in tree hydraulic safety margins predicts the carbon balance of Amazon forests. Nature, 617(7959), 111–117.10.1038/s41586-023-05971-3CrossRef Google Scholar PubMed

Torres-Ruiz, J. M., Cochard, H., Mencuccini, M., Delzon, S., & Badel, E. (2016). Direct observation and modelling of embolism spread between xylem conduits: A case study in Scots pine. Plant, Cell & Environment, 39(12), 2774–2785.10.1111/pce.12840CrossRef Google Scholar PubMed

Tyree, M. T., & Sperry, J. S. (1989). Vulnerability of xylem to cavitation and embolism. Annual Review of Plant Physiology and Plant Molecular Biology, 40(1), 19–36.10.1146/annurev.pp.40.060189.000315CrossRef Google Scholar

US EPA Office of Water (2005). Membrane Filtration Guidance Manual. US Environmental Protection Agency.Google Scholar

Vergeynst, L. L., Dierick, M., Bogaerts, J. A., Cnudde, V., & Steppe, K. (2015). Cavitation: a blessing in disguise? New method to establish vulnerability curves and assess hydraulic capacitance of woody tissues. Tree Physiology, 35(4), 400–409.10.1093/treephys/tpu056CrossRef Google Scholar PubMed

Vespignani, A. (2012). Modelling dynamical processes in complex socio-technical systems. Nature Physics, 8(1), 32–39.10.1038/nphys2160CrossRef Google Scholar

Wason, J., Bouda, M., Lee, E. F., McElrone, A. J., Phillips, R. J., Shackel, K. A. … Brodersen, C. (2021). Xylem network connectivity and embolism spread in grapevine (Vitis vinifera L). Plant Physiology, 186(1), 373–387.10.1093/plphys/kiab045CrossRef Google Scholar PubMed

Watts, D. J., & Strogatz, S. H. (1998). Collective dynamics of ‘small-world’ networks. Nature, 393(6684), 440–442.10.1038/30918CrossRef Google Scholar PubMed

Wheeler, J. K., Sperry, J. S., Hacke, U. G., & Hoang, N. (2005). Inter-vessel pitting and cavitation in woody Rosaceae and other vesselled plants: A basis for a safety versus efficiency trade-off in xylem transport. Plant, Cell & Environment, 28(6), 800–812.10.1111/j.1365-3040.2005.01330.xCrossRef Google Scholar

Yang, D., Pereira, L., Peng, G., Ribeiro, R. V., Kaack, L., Jansen, S. … Tyree, M. T. (2023). A unit pipe pneumatic model to simulate gas kinetics during measurements of embolism in excised angiosperm xylem. Tree Physiology, 43(1), 88–101.10.1093/treephys/tpac105CrossRef Google Scholar PubMed

Zhang, Y., Pereira, L., Kaack, L., Liu, J., & Jansen, S. (2024). Gold perfusion experiments support the multi-layered, mesoporous nature of intervessel pit membranes in angiosperm xylem. New Phytologist, 242(2), 493–506.10.1111/nph.19608CrossRef Google Scholar PubMed

Figure 1. Figure 1 long description.Schematic presentation of the xylem network and pore constriction models. A) The xylem is represented by a 3D cylindrical grid where the row (r$r$), column (z$z$), and angular (ϕ$\phi$) dimensions correspond to the axial, radial, and tangential dimensions of a tree. Each grid cell can either belong to a vessel (red cells) or not (gray cells). Neighboring vessels form intervessel connections in radial direction with probability Pe,rad$\mathcal{P}_{e, rad}$ and in tangential direction with probability Pe,tan$\mathcal{P}_{e, tan}$. B) In each angular slice of the xylem, each cell above a non-vessel cell starts a new vessel with probability NPc$\mathcal{N}\mathcal{P}_c$, while each cell above a vessel cell ends the existing vessel with probability Pc$\mathcal{P}_c$. Note that isolated vessels and dead-ends are removed before further analysis. C) Each pit membrane pore is modeled as a set of constrictions (red), forming a channel through the 3D pit membrane (gray). The smallest constriction diameter Rmin$R_{min}$ defines the effective diameter of each pore, and BPP of a pit membrane depends on the largest Rmin$R_{min}$ among all pores of the membrane.

Table 1. Overview of the abbreviations and symbols used. For values of parameters, see Table 2

Table 2. Xylem network and embolism spreading model parameters for B. pendula

Figure 2. Figure 2 long description.Optimization of probability parameters for xylem network construction. We constructed xylem vessel networks with different combinations of probabilities NPc$\mathcal{N}\mathcal{P}_c$, Pc$\mathcal{P}_c$, Pe,rad$\mathcal{P}_{e, rad}$, and Pe,tan$\mathcal{P}_{e, tan}$ and calculated the vessel tissue fraction (A, B) and grouping index (GI) (C, D) of each simulated network. Two 2D projections of the 4D space spanned by the four probability parameters are shown: for visualizing the effect of NPc$\mathcal{N}\mathcal{P}_c$ and Pc$\mathcal{P}_c$ (A, C), Pe,rad$\mathcal{P}_{e, rad}$ and Pe,tan$\mathcal{P}_{e, tan}$ are fixed to their optimal value, while the effect of Pe,rad$\mathcal{P}_{e, rad}$ and Pe,tan$\mathcal{P}_{e, tan}$ (B, D) is visualized fixing NPc$\mathcal{N}\mathcal{P}_c$ and Pc$\mathcal{P}_c$ to their optimal values. The red square shows the optimal parameter combinations. Note that the z scale of GI (C, D) is logarithmic; the blank cells indicate areas where GI = 0.

Figure 3. SI model replicates the vulnerability curve produced by the physiological model. A) PLC produced by the physiological (red) and SI (gray) models. For calculating the PLC, Keff$K_{eff}$ was averaged across 100 spreading simulations. The transition from normal xylem function to hydraulic failure happened almost step-like over a narrow P$P$ range; for further details, see the main text and Section 4.3. B) The optimized SI spreading probability (β$\beta$).

Figure 4. Figure 4 long description.Behavior of xylem network properties during embolism spreading simulated with the physiological model. Keff$K_{eff}$ and the largest connected component (LCC) size (A, B, C), prevalence (D, E, F), and nonfunctional sap volume (G, H, I) are monitored for P12$P_{12}$ (A, D, G), P50$P_{50}$ (B, E, H), and P88$P_{88}$ (C, F, I) of the physiological spreading model. For each network property, the full line corresponds to the mean across 100 simulations, while the shadowed area shows the standard deviation across iterations. The x axis shows the dimensionless number of simulation time steps that does not directly map to absolute time. Note that the embolism spreading took different number of time steps to saturate at different P$P$ values, and the x axes thus differ between subplots.

Table 3. RMSE between the network property behavior vectors calculated with the physiological and SI spreading models

Figure 5. Figure 5 long description.Behavior of xylem network properties during embolism spreading simulated with the SI model. Similarly as in Figure 4, Keff$K_{eff}$ and the largest connected component (LCC) size (A, B, C), prevalence (D, E, F), and nonfunctional sap volume (G, H, I) are monitored for P12$P_{12}$ (A, D, G), P50$P_{50}$ (B, E, H), and P88$P_{88}$ (C, F, I), the full lines and shadowed areas show mean and standard deviation across simulations, and the x axis shows the number of simulation time steps. Similarly as in Figure 4, note the difference in the x axis maximum values and spacing between subplots.

Figure 6. The SI model produces a vulnerability curve that matches empirical observations. A) The empirical vulnerability curve (red line) of B. pendula is from González-Muñoz et al. (2018) and obtained as a sigmoidal fit to data. The modeled vulnerability curve (gray) is obtained by optimizing the SI spreading probability (β$\beta$) for each P$P$. B) The optimized β$\beta$ values used for constructing the SI vulnerability curve.

Korhonen et al. supplementary material

DOI: https://doi.org/10.1017/nws.2026.10035.sm001

File 626.7 KB

Article contents

Stochastic spreading models reproduce embolism propagation dynamics in angiosperm xylem networks of vessels connected by bordered pits

Abstract

Keywords

Information

1. Introduction

2. Materials and methods

2.1 Xylem network model

2.2 Hydraulic properties of the xylem network

2.3 Embolism spreading

2.3.1 Physiological spreading model

2.3.2 SI model

2.3.3 Vulnerability curves

2.4 Network properties

2.5 Parameter estimation

3. Results

3.1 SI model replicated physiological embolism spreading

3.2 Network properties revealed the embolization process

3.3 SI model reproduced an empirical B. pendula vulnerability curve

4. Discussion

4.1 SI model describes embolism spreading based on a single free parameter

4.2 Xylem network breakdown contributes to hydraulic failure

4.3 Network structure and model parameters affect modeled embolism spreading

4.4 Which compartment model is best suited for modeling embolism spreading?

4.5 Embolism spreading is a directed percolation phenomenon

5. Conclusion

Acknowledgements

Supplementary material

Competing interests

Funding statement

Data availability statement

Ethical standards

Author contributions

References

Korhonen et al. supplementary material

Save article to Kindle

Save article to Dropbox

Save article to Google Drive

Reply to: Submit a response

Your details

You have entered the maximum number of contributors

Conflicting interests