Spatially explicit simulation of deforestation using the Ising-like neutral model

Deforestation is a prominent example of anthropogenic land change. Here we investigate deforestation from a landscape ecology perspective by process-based modeling of a sequence of consecutive landscape patterns starting from a homogeneous forest and ending in homogeneous agricultural land. We propose an Ising-like neutral model to mimic the deforestation dynamics. Our goal is to demonstrate that such a simple conceptualization of the deforestation process is capable of simulating an evolutionary sequence of landscape patterns matching empirical data. The model has only two parameters, one representing an external pressure (for or against deforestation), and another representing a propensity of land for spatial autocorrelation. We perform a suite of simulations to link a process, as parameterized by the values of the two parameters, to the sequence of landscape patterns. We have found a narrow range of parameters for which the model quantitatively and visually matches empirical data. This indicates that the Ising-like neutral model captures the essential aspects of a deforestation process and that the temporal sequence of intermediate landscape patterns does not depend on the details of the process.


Introduction
According to a study by Hansen et al (2013), approximately 2.3 million km 2 of forested land were lost during the period from 2000 to 2012. For regular updates on forest losses, you can refer to the Global Forest Watch online platform at http://globalforestwatch.org. It is important to note that not all of this loss was permanent. Curtis et al (2018) have classified forest losses into the following categories, commodity-driven deforestation, shifting agriculture, forestry, wildfires, and urbanization. Of these five categories, only deforestation and urbanization are permanent. Although, the Food and Agriculture Organization of the United Nations defines deforestation as the conversion of forest to another land use (Shvidenko 2008), in Curtis et al classification the term deforestation is reserved for the conversion of forest to agriculture. Accepting Curtis et al classification, a time evolution of an initially predominantly forested landscape into a predominantly agricultural landscape is the simplest (binary) case of land-use/land cover (LULC) change.
The spatially explicit problem of LULC change can be approached using various methods, as discussed in the review by Ren et al (2019). This review suggests classifying modeling approaches along a spectrum ranging from pattern-to-process methods to process-to-pattern methods. Additionally, these approaches can be categorized based on their methodology, namely empirical, mechanistic, and neutral models (Pearson andGardner 1997, With andKing 1997). Empirical models involve deriving rules for landscape pattern change from historical data and applying them to predict future patterns. Mechanistic models, on the other hand, rely on rules that stem from well-understood and quantifiable biological, economical, or behavioral processes. Lastly, neutral models employ process-neutral rules that generate patterns consistent with observed data.
Most of the literature on deforestation (Hansen et al 2013, Keenan et al 2015, d'Annunzio et al 2015, Feng et al 2016, Li et al 2018, Galiatsatos et al 2020 primarily focus on estimating the magnitude and spatial distribution of deforestation using pattern-to-process empirical methods. In contrast, process-to-pattern methods often take a mechanistic approach, such as economic-based models (Hertel 2018) or agent-based models (Parker et al 2003). However, due to the intricate nature of deforestation, which involves numerous potential drivers and underlying theories (Indarto and Mutaqin 2016), process-to-pattern models are less effective in making accurate predictions compared to empirical models. Moreover, due to their complexity, mechanistic models lack a clear and intuitive connection between the underlying process and the resulting pattern. Such connection can, in principle, be provided by a neutral model if a process-neutral evolutionary rule can be found that is capable of reproducing observed patterns and its parameterization is readily explainable in environmental terms.
In this study, we present and assess a neutral model designed to simulate the evolution of binary landscapes within the context of forest-to-agriculture conversion (FAC). This model is characterized by two key parameters: one representing the strength and direction of external (exogenous) factors, and the other capturing the landscape's propensity for spatial autocorrelation. External factors may include policies related to agricultural land expansion (favoring deforestation) or environmental protection measures (discouraging deforestation). Examples of factors influencing autocorrelation propensity encompass land topography, site accessibility, and zoning. The proposed model offers a notable advantage: it does not require explicit specification of either exogenous or endogenous factors to simulate FAC or comprehend observed dynamics.
The proposed model is a modified version of the Ising model (IM) (Ising 1924, Onsager 1944, Brush 1967, Cipra 1987. Initially introduced as a statistical-dynamical system for understanding magnetic substances (Ising 1924), the IM has become a fundamental concept in statistical mechanics, serving as a paradigm for studying phase transitions across various disciplines beyond physics. It has found applications in diverse fields including social science (Bornholdt andWagner 2002, Stauffer 2008), psychology (Cramer et al 2016, Brandt andSleegers 2021), genetics (Majewski et al 2001), and environmental studies (Ma et al 2019). Within ecology, the IM has been utilized to investigate forest canopy-gap structure (Kubo et al 1996, Kizaki and Katori 1999, Schlicht and Iwasa 2006, model vegetation patterns along regional rainfall gradients in southern Africa (Scanlon et al 2007), and explain emergent phenomena such as masting (Noble et al 2018)-the irregular, periodic synchronous production of large seed crops in perennial plants, where certain years exhibit abundant seed production while others experience limited or no crop.
We define deforestation as a sequential process that begins with a forested site and ends in an agricultural site. When considering sites of substantial size (typically 100 km 2 in our study), this transition unfolds gradually over an extended period, during which the site's landscape evolves from forested to agricultural. To capture this evolution, we rely on a temporal series of site observations, consisting of snapshots that capture the intermediate patterns over time. We refer to this sequential collection of snapshots as an evolutionary track. In our approach, we employ the IM to simulate the site's evolution. The main technical challenge lies in determining the optimal values for the two parameters of the IM model, such that the simulated track aligns closely with the observed track.
Our study aims to address the following key questions: (1) Can the IM be reinterpreted within the broader context of land change modeling, specifically pertaining to the FAC model? (2) What are the various types of evolutionary tracks that can be generated through the application of the IM? (3) How are the values of the IM's parameters, representing the underlying processes, are connected to the outcome (the final pattern)? (4) Is it possible to identify the IM's parameter values that yield a simulated FAC track closely resembling the observed patterns?

Neutral model
Given the distinct nomenclature used in landscape ecology and physics (the origin of the IM), it is essential to establish clear definitions. In this paper, we define a landscape as a rectangular array of cells (small tracts of land), with dimensions L × L. Each cell, denoted as i (i = 1, . . . , L 2 ), is categorized as either forest or agriculture. It is important to note that 'forest' and 'agriculture' refer to specific land cover types, which have broad definitions, particularly in the case of agriculture. The IM, which operates on a two-dimensional square lattice with dimensions L × L, assigns spins to each lattice site. These spins, denoted as i (i = 1, . . . , L 2 ), can take on two values: +1 or −1. Consequently, we can formally transform the landscape into the Ising model by substituting the 'forest' category with +1 and the 'agriculture' category with −1. As the values of cells in the landscape are categorical, they can be represented by various symbols, including numerical values.
When we flatten the 2D array representation of a given landscape, we obtain an equivalent description in the form of a string x 1 , . . . x L 2 with a length of L 2 . This string consists of values of +1 and −1, and in the context of the IM, it is referred to as the system's state or configuration. Therefore, in this paper, the terms landscape, pattern, state, and configuration are used interchangeably, with a preference for pattern or landscape as they are standard terminologies in landscape ecology.
Let ω s = x s,1 , . . . x s,L 2 denote a specific pattern. The IM assigns a numerical value of energy, denoted as E(ω s ), to a pattern or configuration ω s : Here, the symbol ⟨i, j ⟩ indicates the summation over the four nearest neighbors of each cell. In its original application to ferromagnetism, this equation has a solid physical foundation. The more negative the value of E(ω s ), the higher the likelihood of observing the pattern (see Gibbs (1902) or Jeschke (2015) for a contemporary explanation). Thus, the pattern of spins evolves towards configurations with lower values of E(ω s ), while being subject to the constraints imposed by the parameters J and B. In our model, equation (1) is a process-neutral postulate. The values E(ω s ), J, and B are dimensionless and do not possess specific physical meanings, although they can be broadly interpreted. The second part of the postulate states that an initial pattern evolves towards configurations with lower values of E(ω s ). This process is executed using the Metropolis algorithm, which is described in the subsequent subsection. The second term on the right-hand side of equation (1) corresponds to the parameter B, which is interpreted as an extraneous pressure exerted on the landscape. Positive and negative values of B induce cells to change its category from agriculture to forest and from forest to agriculture, respectively. In the first term on the right side of equation (1), the parameter J represents the strength of the land's propensity to maintain spatial autocorrelation. While multiple pairs of parameters (J, B) can lead to an equilibrium state characterized by an agricultural landscape, the evolutionary paths of these states would differ. Our objective is to identify values of (J, B) that produce simulated evolutionary tracks consistent with empirical tracks.

Metropolis dynamics
The IM itself, as described by equation (1), is inherently static. The dynamics is introduced through an algorithm known as the Monte Carlo-based Metropolis algorithm (Metropolis et al 1953). This algorithm iteratively determines the equilibrium configuration for a given system (landscape) with specified (J, B) values. In terms of time, each iteration of the Metropolis algorithm can be considered as a unit, thus establishing the concept of Metropolis dynamics.
The algorithm begins with a given set of parameter values J and B and starts from an initial pattern (a forest). It iterates until an equilibrium pattern corresponding to the specified parameter values is achieved. This process is accomplished through the following steps: (i) Starting from a pattern ω s1 with an energy of E(ω s1 ), a small random change is made to the pattern, resulting in a new, similar pattern ω s2 . (ii) The energy of the new pattern, E(ω s2 ), is computed, typically resulting in a small change from E(ω s1 ).
, the new pattern is accepted since it has lower energy and is thus closer to the equilibrium.
where σ is a random number generated from a uniform distribution between 0 and 1. This criterion allows the algorithm to occasionally accept higher energy patterns to ensure that the global minimum energy pattern is eventually found.
We utilize the (Mathematica 2022) implementation (Nishidate et al 1993) of the Metropolis algorithm for our simulations. Periodic boundary conditions are employed, which means that when a cell along the boundary is randomly selected for a potential change, the neighboring sites are taken from the opposite side of the array. This treatment makes the L × L lattice of cells topologically equivalent to a torus. By using periodic boundary conditions, we mitigate the edge effects in simulations with a relatively small size of L (e.g. 50 cells). To complete one Monte Carlo update, which corresponds to a single time step in the FAC's evolutionary track, the Metropolis algorithm is repeated L 2 times. This ensures that each cell in the array has an equal chance of being selected during a single iteration. Updates are repeated multiple times until the pattern stabilizes, indicating the attainment of equilibrium.

Observations
Our objective is to demonstrate that our model can accurately simulate the evolutionary trajectory of FAC. However, complete direct observations of the entire FAC process at the desired spatial scale are currently unavailable. We find the spatial scale of at least approximately 100 km 2 , referred to as a mesoscale, to be particularly interesting due to the significant environmental impact resulting from deforestation. Therefore, we aim to simulate FAC on a mesoscale to capture the meaningful effects of deforestation on the environment.
The dataset constructed by  encompasses evolutionary trajectories from 1992 to 2015 for over 1.6 million mesoscale landscapes spanning the entire terrestrial landmass . Their analysis revealed that only 22% of these landscapes exhibited noticeable changes and none of them completed the FAC during this period. It follows that the timescale associated with mesoscale landscapes in general, and FAC in particular, extends beyond the available observational record. Hence, we are compelled to rely on observations-based reconstructions.
This situation is not unprecedented and has parallels in other scientific disciplines. For instance, similar challenges have arisen in the study of stellar evolution (Salaris and Cassisi 2005), plate tectonics (Conrad and Lithgow-Bertelloni 2002), and the evolution of life (Koonin 2011), where the timescales of evolution exceed the limits of direct observation for individual objects. To overcome this hurdle, researchers have adopted a strategy of observing numerous objects at different stages of their evolutionary tracks and leveraging these observations to reconstruct complete evolutionary trajectories.
For FAC, such a strategy was used by . They initially classified landscapes in 1992 and 2015 into 64 configurational categories and determined the probabilities of transitioning from each class in 1992 to classes in 2015. By assuming that this probability matrix remains consistent over any 23 year period, they utilized it to reconstruct FAC trajectories. A reconstructed trajectory consists of real landscapes but not from the same site. The outcome of their approach is a collection of reconstructed FAC tracks, sorted in descending order based on their frequency of occurrence. Each reconstructed track represents an ordinal temporal sequence, lacking specific time values between snapshots, ranging from predominantly forested to predominantly agricultural landscapes. We refer to reconstructed tracks as empirical tracks.

Results
Our results can be divided into two parts, (1) connecting the IM's parameters to types of evolutionary tracks, and (2) finding parameters for which IM-simulated evolutionary tracks most resemble stochastic empirical tracks identified by .
The overall design of a single simulation involves beginning the process from an initial condition of an almost homogeneous forest and then evolving the pattern using the IM with specific parameters (J and B). It is important to note that the initial condition is not perfectly homogeneous, as we use a real landscape as the starting point, which contains some non-forest impurities (refer to figures 2 and 3 for visual representations). The landscape is approximated by a square array with a linear size of L = 50 cells. In each iteration of the Metropolis algorithm, a randomly selected cell has the opportunity to change its category, thereby introducing a slight modification to the entire pattern. A total of L 2 = 2500 iterations are executed to ensure that each cell in the landscape has a chance to be modified. After L 2 iterations, the pattern is considered updated and saved for future reference. We conduct up to 500 updates, terminating the simulation when the model reaches equilibrium. Depending on the parameter values chosen, the equilibrium can manifest as a homogeneous forest, transition to a homogeneous agriculture, or evolve into a mixed pattern consisting of both forest and agriculture.
To investigate the relationship between parameter pairs and their corresponding outcomes, we conducted a series of simulations encompassing the following parameter ranges: B ranged from −1.0 to 0.2, while J ranged from 0.1 to 1.0. The range of B incorporates some positive values that promote the preservation of the initial homogeneous forest pattern. However, the majority of values in this range are negative, indicating a preference for the conversion of the initial pattern to agriculture. points where ε ≈ 0 identify simulations that maintain a homogeneous forest pattern. The remaining points identify simulations resulting in a mixed forest/agricultural landscape.
When interpreting figure 1, it is important to keep in mind that B represents the strength of an exogenous force, while J signifies the strength of the landscape's propensity to maintain spatial autocorrelation. Simulations characterized by non-negative values of B (panels 1-3 in figure 1) demonstrate a pro-forest nature, thereby preserving the homogeneous forest landscape unless the value of J is excessively small. In this context, a small J implies a weak autocorrelation propensity, leading to the inability of the landscape to maintain a 100% forest composition despite the external pro-forest pressure.
Simulations characterized by small negative values of B (weakly pro-agriculture) in panels 4-7 of figure 1 result in the formation of the FAC, except when J is either too small or too large. In the regime of small J, where the autocorrelation propensity is weak, forested spots emerge, preventing the landscape from settling into a homogeneous agricultural pattern. The resulting landscape becomes a mixture of agriculture and forest. In the regime of large J, the autocorrelation propensity is so robust that the forest persists despite the external pressure.
Simulations characterized by large negative values of B (strongly pro-agriculture) in panels 8-12 of figure 1 yield the FAC or a near FAC for all values of J within the range of 0.1 to 1. However, it is worth noting that if the value of J is exceptionally large, the forest can still be preserved. Such a regime extends beyond the range of snapshots shown in figure 1 By analyzing figure 1, two additional interesting behaviors emerge. Firstly, for simulations involving negative values of B, a transition occurs from a fully agricultural state to a fully forested state with a small increase in the J parameter. This abrupt alteration in outcome is commonly referred to as a phase transition. Secondly, within the B range of −0.3 to −0.6 (panels 6-9), a hysteresis phenomenon is observed. As J approaches the value of the phase transition from the left, the anticipated response of the landscape shifting from pure agriculture to pure forest lags behind. Consequently, a range of J values exists where simulations can result in either an agricultural or forested outcome, contingent upon the specific evolutionary track. Although hysteresis is also evident for larger B values (panels 10-12) and beyond, it occurs at J values outside the parameter space covered in this study. While the existence of phase transitions and hysteresis aligns with the fundamental principles of the IM upon which our model is based, the detailed investigation of these phenomena falls outside the scope of this paper.

Classification of evolutionary tracks
The results depicted in figure 1 establish a connection between the simulation outcome (the final pattern in a track) and the parameters of the IM (the process). However, our primary focus lies in the nature of intermediate patterns and their dependence on the IM parameters. To investigate this, we conducted calculations for 130 parameter pairs (J, B), following the methodology outlined earlier in this section. In all simulations, an initial landscape representing an almost homogeneous forest served as the initial conditions. Each evolutionary track consists of eight equispaced snapshots: the initial landscape, six intermediate landscapes, and the final outcome landscape. The number of IM updates between each pair of consecutive snapshots remains consistent within a given track but may vary across different tracks, as not all tracks require the same number of updates to attain equilibrium. Figure

Matching simulated tracks to empirical tracks
The upper row of figure 3(A) illustrates the most common empirical FAC track, as identified by . Additionally, another empirical FAC track is shown in the upper row of figure 3(B), which represents a less direct path to FAC, involving several intermediate steps. It is important to note that the time spans between different snapshots are not predetermined and are variable (for more information, refer to ).
In this study, our objective is to identify parameter values for the IM model that result in simulated FAC tracks demonstrating qualitative similarity to the empirical tracks discussed earlier. To accomplish this, we consider pairs of parameters belonging to the Af, Am, and Ac classes (refer to figure 2). Among these pairs, we select the best fit based on the following approach.
Starting from the initial condition represented by the first snapshot in the empirical track, we perform the simulation for a given pair (J, B). After each update, we calculate two metrics characterizing the updated pattern: ε, which indicates the proportion of agricultural cells (composition), and ⟨s i s j ⟩, representing the average correlation between adjacent cells (reflecting landscape autocorrelation and called a texture or a clustering parameter). The Euclidean distance between the metric pairs (ε, ⟨s i s j ⟩) sim and (ε, ⟨s i s j ⟩) emp , corresponding to the simulated pattern and the upcoming snapshot of the empirical track, is then computed.  The updating process halts when the Euclidean distance reaches its minimum, indicating that the simulated pattern has evolved to a configuration most similar to the next snapshot in the empirical track. This pattern is saved as the snapshot in the simulated track, and the procedure repeats until the determination of the last (eighth) simulated snapshot. The best fit pair of parameters is identified by calculating the sum of Euclidean distances between the corresponding simulated and empirical snapshots for all eight time steps. The pair with the smallest sum of distances is considered the best fit, indicating the closest match between the simulated and empirical tracks.
The best fit to the most probable empirical track corresponds to J = 0.7 and B = −0.3, representing the Ac track category. Figure 3(A) provides a comparison of eight snapshots from both the empirical and simulated tracks. It is important to note that the empirical tracks exhibit a stochastic nature, with snapshots serving as archetypes derived from various locations. These snapshots exemplify the expected composition and configuration of intermediate landscapes at different time steps within a specific FAC. However, it should be emphasized that the layouts of snapshots in the empirical track may not be consistent due to the nature of the reconstruction. In contrast, the simulation generates the temporal progression of a single site, resulting in consistent snapshots layouts throughout. To accommodate for the arbitrary orientation of landscapes in the empirical track, certain snapshots within the matching simulated tracks have been translated or rotated on the surface of a torus (as described in section 2) to align with the layout of the landscapes in the empirical track to facilitate visual comparison between the simulated and empirical tracks.
The two tracks exhibit remarkable similarity in terms of their compositional metric (ε) and texture metric ⟨s i s j ⟩. Visually, they also appear highly comparable. It is worth noting that both tracks require more stages (snapshots) to reach a 50% reduction in forest cover compared to the number of snapshots it takes to reach the second 50% reduction. In both the empirical and simulated tracks, the 50/50 split occurs between snapshots five and six. This observation indicates that the simulation accurately captures the acceleration of the deforestation process after the initial half of the forest has already been lost. This acceleration is evident in terms of stages and, as suggested by , it likely corresponds to an acceleration in actual time as well.
The best fit to the second empirical track, which is less common, is J = 0.7 and B = −0.4. It also falls into the Ac track class. It is important to note that this track exhibits several reversals, where deforestation briefly transitions into reforestation (as indicated by a decrease in the value of ε instead of an increase). These reversals may be artifacts associated with the construction of empirical tracks. However, it is worth mentioning that in reality, short episodes of reforestation during the FAC can occasionally occur. To account for these brief reforestation episodes, we employ the model with J = 0.7 and B = 0.4 (a reforestation model) to advance the simulation to the subsequent snapshot. Once again, both quantitatively (in terms of ε and ⟨s i s j ⟩ values) and visually, the empirical and simulated tracks demonstrate a high degree of similarity. The empirical track aligns with the same acceleration pattern observed in the most probable track, where deforestation intensifies after the initial 50% loss of forest cover. The simulation effectively captures this characteristic in the less likely track as well.

Discussion and conclusions
In the Introduction, we presented four key questions that we aimed to address in this study. In this section, we will provide comprehensive answers to these questions based on our results.
To answer the first question, we demonstrated (in section 2) how to construct a neutral model of binary landscape evolution based on the well-known IM. The difference between our model and the original IM is the types of problems they address. The original IM is a prototypical model of phase transitions in statistical mechanics and it is used as such in many fields, including ecology (Pascual andGuichard 2005, Ma et al 2019), to explain observed transitions between organized and disorganized phases. On the other hand, temperature is not a natural variable in LULC landscape evolution, thus, it is not an explicit parameter in our model (formally, we set k B T = 1). A phenomenon of phase transition can still (theoretically) occur, but conditions for its occurrence are indexed by critical values of (J, B) (see section 4) instead of temperature. Overall, we have successfully reinterpreted the original IM to a neutral model of binary landscape's evolution.
The results summarized in figures 1, 2 and S1 provide an answer to the second question. The diversity of tracks that can be simulated by our model is succinctly depicted by the seven different classes of tracks shown in figure 2. Although the initial patterns are the same for all seven categories, subsequent evolution varies depending on the values of the model's parameters J and B. In particular, note that tracks in categories MF and ME lead to a partial loss of forest in the absence of any pressure for expansion of agricultural lands (ME) or even in the presence of small pressure for the preservation of forest (MF). This leads to an interesting theoretical possibility, that under the weak forest conservation policy, some forest can be lost if the propensity for local same-land use aggregation is non-existent or small. However, according 'the first law of geography' (Tobler 1970), naturally occurring landscapes do not have a very fine texture, thus we do not expect such tracks to be observed. Indeed, simulations of the two most frequent empirical tracks belong to the category Ac and are characterized by aggregated patterns (see figure 3).
The third question, concerning the process-pattern connection, is addressed by the classification diagram presented in figure 2. In this diagram, the process is represented by the pairs of parameters (J, B) used in the model. Each class in the diagram groups together pairs of parameters that result in different final equilibrium patterns when starting from the same fully forested initial landscape. This classification diagram establishes the connection between the process (parameter values) and the resulting patterns. Notably, it is possible for multiple classes of the process to lead to the same ultimate outcome, but through distinct evolutionary tracks (as illustrated by the evolutionary tracks under processes Ac and Am in figure 2).
Regarding the fourth question, the results presented in section 4.2 provide an answer. Figure 3 demonstrates our success in identifying values for the parameters J and B that enable the evolution of a forested landscape into an agricultural landscape, with intermediate stages exhibiting patterns that are quantitatively and visually similar to corresponding intermediate stages in two different empirical tracks. Thus, our hypothesis that an Ising-like neutral model can simulate the FAC evolutionary track is supported by our findings. Interestingly, the best fit parameters for the processes explaining the two empirical FAC tracks have values of J and B close to the critical values depicted in panels 6 and 7 of figure 1.
Our overarching hypothesis posits that the Ising-inspired neutral model described in this paper represents the 'universality class' of dynamical processes governing all types of binary LULC landscapes. It is important to note that we employ the term 'universality class' in a broad sense, deviating from its specific definition in physics (Cardy 1996). Nonetheless, this terminology aptly captures our hypothesis, suggesting that binary landscapes exhibit similar evolutionary behaviors, as encapsulated by our neutral model, despite the potential variation in the underlying physical mechanisms across different landscapes.
This paper represents an initial step towards investigating the aforementioned hypothesis, and the obtained results are promising. However, further research is necessary to advance our understanding in this area. This includes conducting additional simulations of FAC and validating the obtained results through observations. In this study, we specifically focused on the evolution of landscapes where the empirical tracks exclusively consist of patterns with forest/agriculture composition at all stages. Moreover, we chose to simulate the entire FAC process, starting from pure forest and concluding with pure agriculture. These choices result in long evolutionary time scales and a limited number of examples. In future studies, we intend to refine our approach by reclassifying LULC maps into binary patterns, distinguishing between focus and background categories. For instance, the forest category can be designated as the focus category category (while also considering other LULC categories as potential focus categories), and the background category can be formed by combining all other categories. Additionally, we plan to explore simulating the FAC process over period shorter than its entire temporal extent. This modified protocol will result in less drastic changes between the initial and final snapshots of evolution. However, such an approach will enable comparisons between simulated tracks and tracks observed in a single site.
Lastly, recall that the IM is specifically designed to simulate the evolution of binary landscapes. However, many real-world landscapes of interest encompass more than two categories. In light of this, it is worth considering the development of a multi-categorical neutral model for landscape evolution in the future. One potential avenue for such development is the utilization of the Potts model (Wu 1982), which serves as a generalization of the IM to encompass more than two components.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary information files).