Estimating process‐based model parameters from species distribution data using the evolutionary algorithm CMA‐ES

Two main types of species distribution models are used to project species range shifts in future climatic conditions: correlative and process‐based models. Although there is some continuity between these two types of models, they are fundamentally different in their hypotheses (statistical relationships vs. mechanistic relationships) and their calibration methods (SDMs tend to be occurrence data driven while PBMs tend to be prior driven). One of the limitations to the use of process‐based models is the difficulty to parameterize them for a large number of species compared to correlative SDMs. We investigated the feasibility of using an evolutionary algorithm (called covariance matrix adaptation evolution strategy, CMA‐ES) to calibrate process‐based models using species distribution data. This method is well established in some fields (robotics, aerodynamics, etc.), but has never been used, to our knowledge, in ecology, despite its ability to deal with very large space dimensions. Using tree species occurrence data across Europe, we adapted the CMA‐ES algorithm to find appropriate values of model parameters. We estimated simultaneously 27–77 parameters of two process‐based models simulating forest tree's ecophysiology for three species with varying range sizes and geographical distributions. CMA‐ES provided parameter estimates leading to better prediction of species distribution than parameter estimates based on expert knowledge. Our results also revealed that some model parameters and processes were strongly dependent, and different parameter combinations could therefore lead to high model accuracy. We conclude that CMA‐ES is an efficient state‐of‐the‐art method to calibrate process‐based models with a large number of parameters using species occurrence data. Inverse modelling using CMA‐ES is a powerful method to calibrate process‐based parameters which can hardly be measured. However, the method does not warranty that parameter estimates are correct because of several sources of bias, similarly to correlative models, and expert knowledge is required to validate results.


| INTRODUC TI ON
The speed and magnitude of projected climate changes are profoundly affecting species distributions, ecological communities and ecosystem processes, and numerous ecological systems are now approaching tipping points (Barnosky et al., 2012;Lenton et al., 2008;Steffen et al., 2018; but see Brook et al., 2013). Large uncertainties on the persistence and the resilience of ecosystems exist. Ecological forecasting has now become a critical tool for managers and decision-makers (Urban, 2015), and robust predictive approaches are necessary to provide reliable projections of species geographic range shifts and ecosystem functioning (Mouquet et al., 2015).
Forecasting the dynamics of ecological systems for the upcoming decades and centuries is very difficult, because ecological systems are extremely complex, influenced by a lot of factors and processes, and climatic conditions with no analogues in the recent past are forecasted to become common (Fitzpatrick et al., 2018;Radeloff et al., 2015;Williams et al., 2007). Ecological models have thus increased in complexity over the last 50 years, incorporating more and more processes described with various degrees of complexity depending on their objectives.
Nowadays, two main types of species distribution models (SDMs) are used to project species range shifts in future climatic conditions: correlative and process-based models (Dormann et al., 2012). The vast majority of currently used SDMs are correlative: they seek to find statistical relationships between various environmental descriptors and species presence and absence. They assume there is an equilibrium between species distribution and environment (equilibrium postulate; Guisan & Thuiller, 2005), that there is no adaptive responses within a generation (no trait plasticity; Berzaghi et al., 2020), and that species niche is stable over macroevolutionary time (niche conservatism; Pearman et al., 2008). Most of them include a fairly large number of predictors (particularly in machine-learning approaches), and consider flexible transformations (linear, quadratic, etc.) and interactions between them (Merow et al., 2014). Even though some authors advocate for 'putting more biology into SDMs' , parameters have no a priori defined ecological meaning (Dormann et al., 2012) and shape of response curves to environmental variables is generally not constrained based on biological considerations. Although these models are not always used correctly (Araújo et al., 2019;Santini et al., 2021), their flexibility makes them an important tool in predictive ecology (Mouquet et al., 2015). They have been widely used especially to generate species range projections under current and future climates (e.g. Guisan & Thuiller, 2005).
Nevertheless, their ability to accurately describe the effects of climate on species distributions has recently been questioned (e.g. Fourcade et al., 2018;Journé et al., 2020;Warren et al., 2021).
For all these reasons, another kind of models has been developed.
Process-based models aim to translate into mathematical equations our knowledge about the physiological and ecological processes involved in an organism's life, such as growth, reproduction, survival, movement and interactions with other livebeings. Process-based models take more time to develop and are more challenging to use, but they might provide a greater comprehension of the complexity of ecosystem dynamics and more robust projections in novel conditions (Evans, 2012;Singer et al., 2016;Urban et al., 2016;Zurell et al., 2016). A wide variety of process-based models exists, from quite simple models (e.g. Kleidon & Mooney, 2000) to much more complex ones (e.g. Dufrêne et al., 2005). They all rely on an explicit representation of processes at stake, with a direct biological interpretation-at least in principle (Connolly et al., 2017). Process-based models may also include some phenomenological relationships on lower-level processes, but never on the pattern itself. The choices about the specific processes to include into the model are made based on theory, empirical observations and the objectives of the research, and modeller subjectivity may play an important role. One of the challenges is to build a model with the appropriate amount of complexity: a too simplistic model might be unrealistic whereas a very complex model could be far beyond our ability to understand it (because of interconnected mechanisms) and calibrate it. Each model relies on different hypotheses with its own balance of complexity, accuracy and parsimony-and thus different numbers of unknown parameters to calibrate. Parameter values of this kind of models are obtained with different methods. Some of these parameters can be individually estimated with field observations or experimental data, or are already available in the literature. When this is not possible, groups of parameters defining a relationship between a process variable and environmental variables or other processes are estimated jointly using inverse modelling methods and data on the processes modelled (e.g. Asse et al., 2020;Cailleret et al., 2020). Calibration (i.e. parameter setting and estimation) is a fundamental step in the modelling process. It allows the model to reproduce the reality with more or less success. The result of the calibration provides insights on the ability of the model to reproduce and explain reality (model predictive power). Calibration of complex models such as process-based SDMs is time-consuming, and modellers are often challenged by the dimension of the parameter space, the complexity of the possible correlations among parameters, and the scarcity of observed experimental data to calibrate them. Therefore, although the philosophy of such models is to measure parameters, statistical inference might be useful when data are not yet available to infer parameter values. Parameter inference can be achieved through several methods which have been developed in the last decades. Most of them fall into two categories: informal and statistical calibration. On the one hand, statistical calibration assumes a data-generating model, and a likelihood function, which quantifies the probability of the observed data given the model parameters, is expressed mathematically.
This likelihood is the foundation of maximum likelihood estimation and Bayesian inference, two commonly used methods for parameter estimation. In practice, deriving the likelihood function can be quite challenging for complex models with many parameters and thresholding. One approach is to use numerical methods such as simulation-based inference or approximate Bayesian computation (e.g. Hartig et al., 2014). However, obtaining a reliable estimate of the likelihood of a complex process-based model can require a large number of simulated samples, which can be computationally expensive. On the other hand, informal calibration uses an informal objective function to measure the discrepancy between the model predictions and the observed data. By minimizing the objective function, one is able to identify the parameter values that best fit the observed data, even though it may not have the same statistical rigour as statistical calibration. This optimization step can be achieved with several methods, either deterministic (e.g. Nelder-Mead method) or stochastic (e.g. simulated annealing, evolutionary algorithms).
Recently, an approach belonging to the evolutionary algorithm family, called Covariance Matrix Adaptation Evolution Strategy (CMA-ES), has been proposed (Hansen & Ostermeier, 2001). One of the advantages of CMA-ES is its ability to cumulate information over iterations in order to adapt its own parameters (in particular the covariance matrix), which makes it more robust to noise.
CMA-ES is especially performant for nonseparable problems (i.e. when the model parameters are dependent) and large search space. This method has been successfully applied in various fields such as aerospace (e.g. Collange et al., 2010), optics (e.g. Gagné et al., 2008) and robotics (e.g. Hill et al., 2020). CMA-ES is acknowledged to be one of the most efficient approaches in continuous black-box optimization  but to our knowledge has never been used in ecology.
Here we explored the feasibility and interests of calibrating process-based SDMs with CMA-ES using species occurrence data as correlative SDMs do (fitted process-based models sensu Dormann et al., 2012). We focused on two forest process-based models of varying levels of complexity to evaluate the ability of CMA-ES to calibrate such models. The two models are PHENOFIT (27-36 parameters; Chuine & Beaubien, 2001) and CASTANEA (77 parameters; Dufrêne et al., 2005). Each model also emphasizes different ecological processes: while PHENOFIT focuses on phenology and how it relates to survival and reproduction, CASTANEA focuses on carbon and water cycles. We focused on three European common tree species, with different range extent and ecological preferences in order to evaluate the algorithms performance in various geographical and climatic conditions. European beech (Fagus sylvatica L.) is one of the most widely distributed broadleaved trees in Europe (from southern Sweden to Sicily and from Spain to northwest Turkey), holm oak (Quercus ilex L.) is an evergreen broadleaved tree native of the Mediterranean region, and silver fir (Abies alba Mill.) is a coniferous tree which mainly occurs in mountain forests of Central Europe and some parts of Southern and Eastern Europe.

| Process-based models
All versions of the models used for this study are coded in Java and distributed by the CAPSIS platform.
PHENOFIT is a process-based species distribution model for forest tree species which focuses on phenology. It relies on the principle that the distribution of a tree species depends mainly on the synchronization of its timing of development to the local climatic conditions (Chuine & Beaubien, 2001). It is composed of several submodels, including phenology models for leaves, flowers and fruits, and stress resistance models. It simulates the fitness (survival and reproductive success) of an average individual using daily meteorological data, soil water holding capacity and species-specific parameters (see Appendix A for details). PHENOFIT has been validated for several North American and European species by comparing their known distribution to the modelled fitness (e.g. Duputié et al., 2015;Gauzere et al., 2020;Morin et al., 2007;Saltré et al., 2013).
CASTANEA is an ecophysiological process-based model which simulates carbon and water fluxes in forests .
The model simulates the ecosystem as an average tree with six compartments (leaves, branches, stem, coarse roots, fine roots and reserves). It is much more complex than PHENOFIT, with several processes described and computed, such as photosynthesis, stomatal opening, maintenance and growth respiration, transpiration and carbon allocation (see Appendix A for details). CASTANEA requires daily meteorological variables and soil characteristics. The model has been initially validated at stand scale for beech , and was then successfully applied to other European species (e.g. Davi et al., 2006;Davi & Cailleret, 2017;Delpierre et al., 2012).
Both models, in their standard version (called here after expert calibration), are parameterized using various sources of information. Some parameters are directly measured or found in the literature such as leaf life span, specific leaf area, LT50 (freezing temperature causing 50% mortality) of leaves, leaf reflectance and so on. Other parameters cannot be measured at all, or their measurement require an enormous effort that cannot be deployed for a large number of species. These parameters are thus inferred by inverse modelling using either Bayesian methods or optimizing methods and data on the specific process they are describing. This is for example the case for phenology models for which all parameters cannot be measured, especially those describing bud dormancy break regulation, since no method allow to measure dormancy break precisely so far (except for a few fruit tree species). Finally, a few parameters are prescribed based on expert knowledge as no data to estimate them exist. For this study, the standard versions of both models were run for the three species using the same climatic data used to do the inverse calibration (see Section 2.2.1), in order to compare the correctness of the predictions obtained with the two types of calibration.

| Climate and soil data
Raw climatic variables were extracted from ERA5-Land hourly dataset (Muñoz Sabater, 2019, 2021) from 1970 to 2000, at a spatial resolution of 0.1 degree in latitude and longitude. We calculated the daily mean values of the following variables used by PHENOFIT and CASTANEA: minimum, mean and maximum daily temperatures, mean dewpoint temperature, daily precipitation, daily global radiation and daily mean wind speed. We computed the daily relative humidity with the ratio of vapour pressure and saturation vapour pressure (both calculated with Clausius-Clapeyron equation) using humidity r package (Cai, 2019). Daily potential evapotranspiration was calculated with Penman-Monteith equation (FAO standard of hypothetical grass reference surface) using a slightly modified version of the ET() function in Evapotranspiration r package (Guo et al., 2016).
Water content at field capacity and wilting point data were extracted from EU-SoilHydroGrids (Tóth et al., 2017) which is at 1 km resolution. Percentage of sand, silt and clay particles, percentage of coarse fragments, bulk density and soil depth were extracted from SoilGrids250m (Hengl et al., 2017) at a 250 m resolution. These data (except for soil depth) are provided at seven soil depths, so we summarized them (weighted sum or weighted mean) taking into account each layer width and total soil depth. Finally, all variables were upscaled at the ERA5-Land spatial resolution 0.1 using bilinear interpolation.
2.2.2 | Tree occurrence data used for the calibration Sources of occurrence data are known to differ even for common European trees (Duputié et al., 2014) and this makes it quite challenging to gather comprehensive data at a sufficient spatial resolution all over Europe. The occurrence data we used essentially rely on the EU-Forest dataset (Mauri et al., 2017) which benefits from inventory and monitoring programmes implemented in most European countries. As EU-Forest is limited to forest ecosystems, we combined it with presence records extracted from the Global Biodiversity Information Facility (GBIF, 2022, see Appendix B for all download links) but removing observations outside natural species ranges as defined by Atlas Flora Europeae (AFE ;Jalas & Suominen, 1972 and EuroVegMap (Bohn et al., 2003). By doing so, we also included occurrences of isolated native trees living outside forests, excluding records from arboreta or gardens where the species would have been planted as an exotic. For holm oak, we also added occurrence records in the Mediterranean Basin from the WOODIV database (Monnet et al., 2021), leaving out EU-Forest and GBIF records we had already gathered. We upscaled all species records at the ERA5-Land resolution (0.1, see Section 2.2.1), that is, the species is considered to be present in the cell if there is at least one record. We finally obtained 21,458 occurrence cells for beech, 6653 for holm oak and 5385 for silver fir (see Appendix B for details).
All the datasets described above are presence-only data.
Therefore, we generated cells where species are supposed to be absent, that is, pseudo-absence cells. In order to avoid as far as possible creating false absence data, we used EU-Forest cells where the species is not reported present as pseudo-absence cells. We assumed that national forest inventories were exhaustive (which is not true since only specific forest plots in a 0.1° cell are monitored). We obtained 25,423 absence cells for beech, 37,931 for holm oak and 38,365 for silver fir (see Appendix C).
We selected subsets of 2000 points (1000 presences and 1000 pseudo-absences) in order to reduce computational costs. For each species, we generated 10 presence clusters (k-means algorithm) of similar bioclimatic conditions based on annual climate normals computed with r package dismo (Hijmans et al., 2021) and ERA5-Land variables. In each cluster, we randomly sampled a number of cells where the species is present proportional to the total number of a number of cells where the species is present in the cluster. The aim of this stratified random sampling was to make sure that all species environmental preferences were proportionally represented. We then randomly sampled the same number of pseudo-absence cells (see Appendix B for details).
Regarding PHENOFIT model, we calibrated 10 times each species parameter set, with five repetitions on two random subsets of presences/pseudo-absences, except for beech. In the latter case, we ran 10 repetitions on 10 subsets (i.e. 100 calibrations) to investigate both the effect of subsampling and the effect of stochasticity.
Since CASTANEA computing time was much higher (see Table 1), we ran only two calibrations for each species (on two different random subsets).  (Hansen, 2006;Hansen & Ostermeier, 1996;Hansen & Ostermeier, 2001). It is based on the principle of evolutionary biology, via recombination, mutation and selection of the most fit candidate solutions (i.e. parameter sets providing the best predictions).

|
At each iteration: • candidate solutions are evaluated, that is, model runs times with different parameter sets and the objective function is evaluated.
• the best candidate solutions are selected.
• the weighted mean candidate solution m is computed (mean of the best parameter sets weighted by their objective function value).
• covariance matrix C and step size are updated (with information accumulated over several consecutive iterations).

| CMA-ES in practice
One of the advantages of CMA-ES is that it does not require a complex parameter tuning: as best parameter values at a given time of the optimization process might no longer be efficient later,

CMA-ES implements an internal adaptation of its parameters.
We only chose the number of candidate solutions , depending on the optimization problem complexity ( was set to ∕ 2). The default recommended value for is 4 + 3ln(N), where N is the number of parameters to calibrate (i.e. ∈ 14, 17 in our case). We set = 20, in order to improve the global search capability (Hansen & Kern, 2004) and take advantage of the computation power at our disposal. All model parameters were linear scaled into 0; 10 so that the same standard deviation can be applied to all parameters: here we chose = 2 (see Nikol aus Hanse n perso nal website for practical hints on variable encoding). Our stopping criterion for the optimization procedure was the budget, that is, the number of model runs.

| Calibration results
Calibrations using species distribution data are hereafter called inverse calibrations, and calibrations based on expert knowledge, observations and measurements of the processes modelled are called expert calibrations.

CMA-ES calibration of PHENOFIT model allows an average
17.2% increase in AUC across the three species compared to expert calibration ( Figure 1). The maximum increase is obtained for silver fir, from 0.72 to 0.9 (25%).
CMA-ES calibration of CASTANEA allows an average 23.7% increase in AUC compared to expert calibration (Figure 2), and a maximum increase obtained for holm oak (34.7%).

| Variability of calibration performance
The 100  no subset led to an overall better prediction of the species distribution (see Figure 3b). The total AUC ranged from 0.881 to 0.914, with a mean value of 0.896.

F I G U R E 1
Species distribution maps obtained with PHENOFIT expert and inverse calibrations, compared with observed species occurrences. Optimal threshold to dichotomize model predicted fitness index in presence/absence is the Youden index-based cut-off point. Note that models predict species climatic niche which is larger than the realized niche that corresponds to species presence map.

| DISCUSS ION
Our results showed that CMA-ES is an efficient optimizer for inverse calibration of complex ecological models. The algorithm was able F I G U R E 2 Species distribution maps obtained with CASTANEA expert and inverse calibrations, compared with observed species occurrences. Optimal threshold to dichotomize model predicted carbon reserves in presence/absence is the Youden index-based cut-off point. Note that models predict species climatic niche which is larger than the realized niche that corresponds to species presence map. Note also that CASTANEA cannot be used in high-latitude regions (grey area).
to find parameter sets that significantly improved the predictions of the calibrated models compared to the expert parametrization.
However, our study also highlighted the issue of nonidentifiability of parameter values due to the data limitation and the dependence between process-based model parameters, which may result in diverging parameter values even though the calibrated models describe the observed species distribution well.

| Performance and advantages of CMA-ES to calibrate complex models in ecology
Here we demonstrated that inverse calibration with CMA-ES is feasible and provide good results for complex and runtime-expensive ecological models.
With a subsampling of species occurrence data, the algorithm succeeds in finding parameter sets which provide higher AUC values. The predictions of the calibrated models are sharply improved compared to expert parametrization (Figures 1 and 2 CMA-ES is a 'generic' optimizer which can be applied to various problems. It is easy to use as it does not require an extensive tuning to efficiently explore the parameter space. We only had to choose the number of candidate solutions , and the initial search region (initial starting point and step size ). As well as being quasi parameter-free, CMA-ES has several structural advantages particularly useful for complex optimization problems. First, the algorithm's covariance matrix enables the learning of second-order information, which provides insights into pairwise dependencies between parameters. Second, the covariance matrix adaptation of CMA-ES is highly efficient in handling ill-conditioned and nonseparable problems (Hansen et al., 2011). Last, CMA-ES update mechanism of step size (i.e. the mutation force) helps prevent premature convergence (Hansen & Ostermeier, 2001), allowing the algorithm to explore more of the search space. CMA-ES has been shown to outperform several other optimization algorithms , and is usually the most efficient method when the target cost (i.e. the number of objective function evaluations) is about 100 × N (N being the dimension of the parameter search space, Bäck et al., 2013).

F I G U R E 3
Effects of data subsampling and stochasticity on CMA-ES calibration using the PHENOFIT model for beech: (a) calibration AUC (calculated only with calibration cells) and (b) total AUC (calculated with all presence/absence cells). Each colour is a different subsampling of occurrence data, each point is a calibration run. Diamonds (with black border) are mean AUC values. On (a), the grouping letters represent the multiple comparisons with pairwise Dunn's tests.

| Nonidentifiability of parameter values
Equifinality is a common problem encountered during model calibration, where multiple parameter sets can produce equally good fits to the observed data. This issue can arise due to several factors, including limitations in the quantity or quality of available data, competing processes within the model, and parameter interactions that affect the model output in complex ways.
In both models used in this study, biological mechanisms are explicitly calculated in several submodels (e.g. a leaf unfolding submodel or a stomatal opening submodel). A submodel output has inevitably a significant influence on the other submodels as biological processes can be highly dependent with feedbacks: in CASTANEA, for example, the stomatal opening affects the photosynthesis, and

| Methodological issues and perspectives
Our goal here was to investigate the performance of CMA-ES to calibrate quite complex process-based species distribution models using species occurrence data. We did not attempt to validate our parametrizations using temporally or spatially independent data, In addition, in our case here, land use management probably also play an important role in shaping tree realized distribution, while not being addressed in the models. We included GBIF occurrence data, and especially as much as possible isolated native tree records outside forests, to help correct this problem, but it is impossible to be exhaustive at the spatial resolution of 0.1°. At this scale, local variations of soil characteristics and of competitive interactions among trees (e.g. along an altitudinal gradient) can also not be considered.
Finally, several authors advocate for process-based modelling approaches relying upon species response functions that are a priori defined (e.g. Higgins et al., 2020). However, the main limitation of such models is the data availability to infer their parameters (Urban et al., 2016). Expert parameterization is often long and arduous. One possible way to facilitate parameter value estimation would be to use inverse calibration, and we demonstrated here that CMA-ES can be a powerful optimizer to this end. For example, CMA-ES driven by species occurrence data could be used to calibrate submodels whose parameter values cannot be measured or are too hard to measure experimentally.
However, when a structural correlation exists (i.e. trade-off between processes that are inherently present and interconnected in the model, as in the leaf unfolding submodel), inverse calibration might not provide the right parameter estimates. In such a case, expert knowledge, observations and measurements are necessary to determine a posteriori which estimates are the most realistic. This is possible in the case of process-based SDMs, and usually not feasible in the case of correlative SDMs. A combination of both expert and inverse calibrations might offer a new perspective for spreading the use of process-based models in predictive ecology, especially for climate change impact studies.

Isabelle Chuine devised the main conceptual ideas. Victor Van der
Meersch worked out the technical details, performed the numerical calculations and wrote the first draft of the manuscript. The two authors discussed the analyses and the results, and contributed to the final manuscript.

ACK N O WLE D G E M ENTS
The authors thank Hendrik Davi for helping us in using the CASTANEA model. We are also deeply grateful for many helpful comments from Florence Tauc. We also thank François de Coligny, manager of the CAPSIS platform, Gilles Le Moguedec and the GenOuest and TGCC teams for their support. Finally, we thank Florian Hartig and another anonymous reviewer whose comments and suggestions helped us improve and clarify this manuscript. Victor Van der Meersch was supported by a GAIA doctoral school PhD Fellowship.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors have no conflict of interest to declare.

PE E R R E V I E W
The peer review history for this article is available at https: