# A workflow to discover partial differential equations from data: Application to the dynamics of tree biomass

Emilie Peynaud, Paulin Melatagia, Serge Stinckwich, Jean-François Barczi

► To cite this version:

Emilie Peynaud, Paulin Melatagia, Serge Stinckwich, Jean-François Barczi. A workflow to discover partial differential equations from data: Application to the dynamics of tree biomass. MethodsX, 2025, 15, pp.103560. 10.1016/j.mex.2025.103560. hal-05287421

HAL Id: hal-05287421

https://hal.science/hal-05287421v1

Submitted on 29 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A workflow to discover partial differential equations from data: Application to the dynamics of tree biomass

Emilie Peynaud<sup>a,b,*</sup>, Paulin Melatagia<sup>c</sup>, Serge Stinckwich<sup>d</sup>,
Jean-François Barczi<sup>a,b</sup>

## GRAPHICAL ABSTRACT

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

**Keywords:**

Physics informed neural network
Parameter estimation
Data-driven modelling of plants
Theory-guided data model

## ABSTRACT

Mixed data and theory driven methods are promising approaches that can be used to bring better understanding of complex dynamics in life sciences. For vegetation growth, integrated knowledge may be lacking to design theoretical models like partial differential equations (PDE). This lack can be complemented by using data. The method presented in this paper is a generic computational workflow called CEDI that aims at discovering PDE models from data. As an illustration, we tested the workflow on biomass dynamics of three different 3D trees of specific architectural types.

<footnote><sup>a</sup> CIRAD, UMR AMAP, F-34398 Montpellier, France</footnote>

<footnote><sup>b</sup> AMAP, Univ Montpellier, CIRAD, CNRS, INRAE, IRD, Montpellier, France</footnote>

<footnote><sup>c</sup> Université de Yaoundé I, Cameroon, IRD, UMMISCO, F-93143, Bondy, France</footnote>

<footnote><sup>d</sup> United Nations University Institute in Macau, Macau SAR, China</footnote>

<footnote>Related research article: None. For a published article: None.</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: emilie.peynaud@cirad.fr (E. Peynaud).</footnote>
* The name CEDI represents the four steps composing the workflow: data Collection, Extrapolation, Differentiation and Identification.

* The originality of this workflow is twofold: first, it encompasses the whole modeling process from the definition of the variables to the design of a PDE, and second it has been designed to be generic in a sense that it can apply to any dynamics and it covers most existing data driven PDE discovering methods.

* The workflow offers a framework to better understand data driven PDE discovering methods and a tool for modeling any dynamics, provided that right data and knowledge and also good algorithm settings are available.

## Specifications table

<table><tr><td>Subject area</td><td>Bioinformatics</td></tr><tr><td>More specific subject area</td><td>Partial differential equations and neural networks</td></tr><tr><td>Name of your method</td><td>Workflow for data Collection, Extrapolation, Differentiation and Identification for partial differential equation discovery</td></tr><tr><td>Name and reference of original method</td><td>Physic Informed neural network</td></tr><tr><td>Resource availability</td><td>Dataset: J.-F. Barczi [1], “Plant growth simulation data”, https://doi.org/10.18167/DVN1/VK0J16, CIRAD Dataverse, V1.<br>Source code: https://gitlab.cirad.fr/amac/personal/emiliepeynaud/pydata2pde.git<br>Software Hash IDentifier: swh:1:dir:fd7172e67af912dddd7f2f50cc6f42f7adb913f</td></tr></table>

## Background

Modeling approaches are crucial tools to characterize tree growth and dynamics in complex plant systems, such as cocoa-based agroforestry systems. Partial differential equations (PDE) are relevant models to represent space and time variations that govern the dynamics of any phenomenon, including those in life sciences. PDEs benefit from a well-established theory that mathematically legitimates this formalism. It applies to various space and time scales, which is a real advantage when dealing with plant and vegetation dynamics. However, estimation of parameters of PDE models still raises questions. The purpose of the workflow presented here is to characterize any dynamics by discovering the PDE that governs the underlying phenomenon using data. The conventional way to derive a PDE model is to formulate with mathematical expressions the known principles that govern the dynamics, considering a given set of assumptions. The workflow detailed in this paper gives an alternative within the paradigm of theory-guided data science [13] that uses data to compensate for a potential lack of knowledge on the principles governing trees' dynamics.

Methods for discovering PDE models from data are gaining a growing interest especially in computational physics [7]. From measurements of the PDE variable at some points in space and time, these methods consist in computing additional estimates of the variable and its partial derivatives to build a library of candidate terms. Among this library, the methods then select a parsimonious combination of terms that best describes the dynamics of the variable. This identification step can be seen as a particular optimization problem which is solved thanks to techniques like statistical inferences [10], sparse regressions [17], sparse least square methods [18] or physic informed neural networks [3,16,19,6]. Most of the work cited here is focused on the optimization problem and is based on synthetic data computed from direct solving of a known PDE and it deals with application in physics where the constitutive laws and hence the PDE forms are well established [14]. This literature review showed that all cases shared common methods for addressing this PDE identification problem that may be encapsulated into a generic workflow providing guidelines for similar cases.

The way to discover PDE proposed in this paper takes the form of a workflow called CEDI and is designed to be generic and adjustable for applications where knowledge or data is lacking. Its originality lies in the embedding of the whole modeling process which is splitted into four formal steps: 1) the definition of the nature of the variable and the data collection (C), 2) the estimates of the values of the variable (E like extrapolation) and 3) its derivatives (D) and 4) the identification (I) of the PDE terms. This gives to the workflow its modularity and its genericity which are the main changes over existing works. For the sake of illustration and as a proof of concept, we applied the workflow to the discovery of PDE governing the biomass dynamics of growing trees from realistic data, using neural networks in the extrapolation and identification steps. The biomass dynamics are considered here as the variations of the leaf and wood mass over space and time. The application to trees as individuals or as communities is not straightforward, but our workflow can complement existing experimental and modeling approaches to better understand the dynamics of plants. The workflow may open new perspectives to characterize dynamics of plants through data by selecting the appropriate differential operators according to both knowledge and measurement of the dynamics. It may also help to give biological meaning to differential operators versus plants by performing successive runs on contrasted dataset.
Method details

This section consists of two parts. The first one is an overview that describes the generic CEDI workflow, which is the main topic of this paper. The second one consists of a case study to apply the CEDI workflow to the specific domain of tree growth PDE modeling, providing a practical implementation to demonstrate and clarify the proposed workflow.

Overview of the CEDI workflow and genericity

The purpose of the generic workflow presented in this paper is to characterize any dynamics by discovering the PDE that governs
the underlying phenomenon using data. By dynamics, we mean the way in which the system under study evolves over time and space
under the influence of the processes to which it is subjected. Let us assume that the phenomenon of interest can be monitored or
described through the measurement of a variable that varies with respect to time and space. We assume that this variable has sufficient
regularity (such as continuity and differentiability) and satisfies a PDE defined in a bounded domain, together with appropriate initial
and boundary conditions. The PDE also involves parameters that quantify the causality of the dynamics. The workflow consists in four
main sequential actions performed on data that migrates from one step to the next: the data collection (C), the extrapolation (E), the
computing of the partial derivatives (D) and the identification (I) of the PDE terms as illustrated in Fig. 1. The purposes of each of the
four formal steps are defined below:

Data collection

The first step is to gather the knowledge available about the phenomenon according to the scientific question motivating the modeling process. Then one determines the traits that vary over space and time and that are meaningful with respect to the scientific question. Let us consider only one trait to simplify the presentation. This trait is mathematically represented by the variable of the PDE and its values must be amenable at various space and time points to build a dataset that faithfully shows the dynamics, as mentioned in Fig. 1. The construction of the dataset depends on the scientific fields and on the techniques available. For example, it may be built thanks to direct measurements of the trait or thanks to the reconstruction of chrono-sequences from a single observation. In this paper, we apply the workflow on measurements of biomass on simulated growing trees across time and space. The dataset may also be built thanks to automatic tools that extract a meaningful dataset from a large and heterogeneous database. Temporal and spatial sampling are supposed to guarantee the representativeness of the dynamics in a dataset of reasonable size and for a reasonable cost. How to know if the dataset faithfully describes the dynamics is a key question that has no general answer. For each case, it may be answered based on requests of both experts on the studied subject and experts on numerical processing.

Extrapolation

This step consists in the design of a surrogate model able to mimic the dynamics of the variable of interest to smooth or complete the
data previously collected. Thanks to the surrogate model, one can compute estimates of the variable at any space-time point for
example in cases of geographically missing data (if the sampling was not exhaustive enough across the experimental domain) or time
missing data (if the sampling time step was not fine enough). This step requires the implementation of an interpolation method
 polynomials, splines, neural networks) on which the next step is based for the computation of missing values of the variable but also
estimates of the partial derivatives.

Differentiation or computation of the derivatives

Thanks to the surrogate model, the third step consists in computing estimated values of the partial derivatives at some space-time
points, not necessarily those of the collected data, in order to build the library of candidate terms for the PDE (also called catalogue in
the following). The candidate terms are chosen according to the available knowledge about dynamics. This selection may be based on

<figure><img src="image_4.png" /><figcaption>Fig. 1. The CEDI workflow for PDE discovery from data.</figcaption></figure>
domain knowledge or on empirical observations. The type of surrogate model often determines the method used to compute the derivatives.

## Identification

The goal of the last step is to select the candidate terms of the PDE that best represent the dynamics and their degree of contribution to it. This is formulated as an optimization problem that minimizes a cost function defined as the residual of the PDE. Knowledge about dynamics helps to design the a priori form of the PDE made of a weighted sum of partial derivatives. If the set of partial derivatives involved in the PDE is known in advance, the identification step reduces to a parameter estimation problem where the goal is to estimate the weight of each partial derivative. If there is a lack of knowledge about dynamics, one can consider several partial derivatives as candidates (they will form a library of partial derivatives or catalogue). Then, selection methods allow us to select the limited number of terms of the PDE. The identification step requires special attention since it raises mathematical questions (well-posedness) as well as implementation and numerical questions (convergence), but also modeling questions (representativeness of the resulting PDE).

The main strength of the CEDI workflow lies in its generalizability in the sense that it encompasses every study about data-driven PDE extraction methods that we found in literature. This is illustrated in Table 1, where we have split some examples of these methods to fit the four steps of the CEDI workflow. These examples let us suppose that CEDI workflow will help any team to prevent and overpass standard issues that rise while discovering PDE from a dataset.

## Detailed presentation with application to tree biomass dynamics

In the following we detail the application of the workflow to the modeling of tree biomass dynamics for the sake of illustration of PDE discovering process from data. For this purpose we defined a very simple Multi Layer Perceptron (MLP) architecture strong enough to provide promising results. Further studies on the particular plant growth question may need more accurate neural networks for step E and I of CEDI and an improved dataset.

## Data collection

Usually, the modeling process starts with the definition of a scientific question and state of the art on the subject. Here the goal is to illustrate the CEDI workflow, so we intentionally simplify what a tree is and how it behaves through space and time. Let us assume the following. The dynamics of tree growth can be seen as the result of continuous biomass creation through photosynthesis (occurring in

<table><caption>**Table 1**
Some examples of data-driven PDE extraction methods decomposed into the CEDI workflow.</caption><thead><tr><td></td><td>Data Collection (C)</td><td>Extrapolation (E)</td><td>Differentiation (D)</td><td>Identification (I)</td></tr></thead><tbody><tr><td>[11]</td><td>Data computed on a grid by adding a Gaussian noise to the analytic solution of an advection equation. Real data application to the prices of SP-500 European call options</td><td>B-splines tensor product approximation introduced as a penalization term in the identification step</td><td>Differentiation based on the B-spline decomposition through the B-spline derivatives</td><td>Frequentist inference based on the maximisation of a penalized least square criterion and Bayesian inference where the frequentist least square criterion is translated into a prior distribution</td></tr><tr><td>[18]</td><td>Data sampled and noised on a grid from simulation performed by direct solving of the PDE on a higher resolution grid</td><td>Fourier transform</td><td>Spectral method based on the Fourier transform for the space derivatives, and first-order backward difference method for time derivatives</td><td>L1-regularized (sparse) least-square optimization problem solved with the Douglas-Rachford algorithm</td></tr><tr><td>[17]</td><td>Random sampling and noising of spatial points and uniform sampling for time points from simulated data by direct solving of the PDE</td><td>Polynomial interpolation</td><td>Finite difference technique applied to a localized points near each measurement points and based on the polynomial interpolation</td><td>Sequential threshold ridge regression (STRidge) consisting in a ridge regression problem with hard threshold</td></tr><tr><td>[16]</td><td>Random sub-sampling and noising of high resolution dataset generated by direct solving of the PDE</td><td>Deep neural network</td><td>Automatic differentiation (TensorFlow)</td><td>Training of the deep neural network (of the extrapolation step) with the L-BFGS algorithm and residual of the a priori known PDE included in the loss function</td></tr><tr><td>[3]</td><td>Measurement of temperature distribution from weather stations in Sweden. Data augmentation by linear interpolation in space and time. Scaling of the data.</td><td>Deep neural network trained with BFGS or L-BFGS algorithms</td><td>Back-propagation or automatic differentiation (TensorFlow or PyTorch)</td><td>Feedforward neural network trained with gradient based methods with L1 regularization for sparsity and feature selection techniques</td></tr><tr><td>[19]</td><td>Generated by direct solving of the PDE with Chebfun package in MATLAB. 1D in space with noise and limited data</td><td>Rational Neural Network (RatNN) with Adam and LBFGS optimizers</td><td>Autodiff capabilities provided by PyTorch</td><td>1) Rational Neural Network with Adam and LBFGS optimizers and re-sampling of the collocation points. 2) Sparse regression technique based on recursive feature elimination (RFE) algorithm</td></tr></tbody></table>
the leaves) and biomass allocation (through phloem flux). The tree biomass is distributed among the plant organs: the leaves, the branches, the roots and the fruits. All along its life, a tree never stops creating and losing new organs so that its biomass is constantly changing over time and space. These dynamics are not well known at the scale of the plant nor the plot, but the above-ground biomass of a tree can be non-destructively monitored. The purpose here is to see how the CEDI workflow could help to characterize the dynamics of the above-ground tree biomass. So, from this angle of view, we define the variable of the PDE as the above-ground biomass of a tree and we note it b, like biomass.

Mathematically speaking, the variable b varies over a space domain $\Omega$ in $\mathbb{R}^3$ encompassing the tree and over a time frame $[0, T]$ with $T > 0$, where $T$ is the lifespan of the tree. That is, we consider $b$ as a function $b: \Omega \times [0, T] \to \mathbb{R}$. We assume that $b$ satisfies sufficient regularity assumptions. We denote $x = (x_1, x_2, x_3) \in \mathbb{R}^3$ the space coordinates in $\Omega$ associated with the Cartesian frame $\{\vec{e}_1, \vec{e}_2, \vec{e}_3\}$ in $\mathbb{R}^3$ and $t$ the time coordinate between 0 and $T$. Assuming that $b$ is differentiable with respect to $x_1, x_2$ and $x_3$ (up to order $k$) and denoting $\gamma = (\gamma_1, \gamma_2, \gamma_3) \in \mathbb{N}^{+3}$ a multi-index of order $k = |\gamma| = \gamma_1 + \gamma_2 + \gamma_3 \ge 1$, we introduce the notation $D^r b(x) := \frac{\partial^r b(x)}{\partial x_1^r \partial x_2^r \partial x_3^r}$ of a $k^{th}$-order partial derivative of $b$ and $D^k b(x) := \{D^r b(x) \text{ such that } |\gamma| = k\}$ the set of all the $k^{th}$-order partial derivatives of $b$. We assume that the dynamics of $b$ can be described by a partial differential equation of the form:

$$ \frac{\partial b(x, t)}{\partial t} = F_{\beta}(D^k b(x, t), D^{k-1} b(x, t), \dots, Db(x, t), b(x, t), x), \quad \forall x \in \Omega \text{ and } 0 < t \le T, \qquad (1) $$

where the right hand side is an application $F_{\beta}: \mathbb{R}^{3k} \times \mathbb{R}^{3k-1} \times \dots \times \mathbb{R}^3 \times \mathbb{R} \times \Omega \to \mathbb{R}$ that depends on $n$ real parameters stored as a vector $\beta \in \mathbb{R}^n$ with $n \in \mathbb{N}^+$ and $n \ge 1$.

Denoting by $\partial\Omega$ the boundary of $\Omega$, the function $b$ also satisfies boundary and initial conditions respectively given by $b(x, t) = 0, \forall (x, t) \in \partial\Omega \times ]0, T]$ and $b(x, 0) = b_0(x), \forall (x, t) \in \Omega \times \{0\}$ where $b_0$ is a known function defined from $\Omega$ to $\mathbb{R}$.

At this stage of the workflow, we do not know the true form of the right-hand side of Eq. (1). If the list of partial derivatives is a priori known, the problem reduces to calibrate the PDE model, that is to find the value of $\beta$. According to [5] and [4], possible candidate terms for $F_{\beta}$ could be the biomass and its first and second order partial derivatives which gives the classical advection-reaction-diffusion operator:

$$ \frac{\partial b}{\partial t} = \rho b - v \cdot Db + D \cdot (Db) \text{ in } \Omega \times ]0, T] $$

With $D := \left( \frac{\partial}{\partial x_1}, \frac{\partial}{\partial x_2}, \frac{\partial}{\partial x_3} \right)^T$ the gradient vector, $\rho \in \mathbb{R}$ the parameter of reaction (or growth rate), $v = (v_1, v_2, v_3)^T \in \mathbb{R}^3$ the advection (or transport) vector which is homogeneous with a velocity and $d = \text{diag}(d_1, d_2, d_3) \in \mathcal{M}^{3 \times 3}(\mathbb{R})$ the matrix of diffusion where $d_i \in \mathbb{R}^+$ for $i = 1, 2, 3$. In that case, there are seven real valued parameters $\beta = (\rho, v_1, v_2, v_3, d_1, d_2, d_3)^T$ to determine. In the general context of tree biomass, we do not have access to these parameters' values through direct measurements, but we do have estimates of the biomass $b$ that can be used to deduce the values of the parameters $\beta$.

Still for the sake of illustration and as a proof of concept of the workflow, we choose to build a biomass dataset from 3D digital mock-ups simulated by the dedicated software AmapSim [2]. This software aims at representing plant growth according to botanical knowledge on the fine architecture of the trees [8] and its variability (in the sense of [12]). The software generates 3D plant shapes (Fig. 2, Fig. 3 and Fig. 4) accurate to reality in terms of organ numbers and organization over time. Based on those mock-ups, we virtually performed a set of $p \in \mathbb{N}^+$ measurements of the biomass. Each measurement consists in a triplet $(b_i^m, x_i^m, t_i^m)$ where $b_i^m$ is the $i^{th}$ biomass value at the $x_i^m$- coordinates at time $t_i^m$. The resulting dataset is available at [1].

<figure><img src="image_5.png" /><figcaption>Fig. 2. Tree growth of Massart type. Mock-ups obtained by the AmapSim software.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 3. Tree growth of Prevost type. Mock-ups obtained with the AmapSim software.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Tree growth of Rauh type. Mock-ups simulated by the AmapSim software.</figcaption></figure>

### Extrapolation

In view of the computation of the parameters $\beta$, we need a tool, let say a surrogate model, to compute extrapolated values of the biomass at any space-time points. Let $b_{\alpha}^s$ denotes that surrogate model which depends on parameters $\alpha$. The surrogate model $b_{\alpha}^s$ is built such that its parameters $\alpha = (\alpha_1, ..., \alpha_N) \in \mathbb{R}^N$ minimize the following loss function

$$J_1(\alpha) = \frac{1}{2} \sum_{i=1}^{p} |b_{\alpha}^s(x_i^m, t_i^m) - b_i^m|^2 + \frac{c_{\alpha}}{2} \|\alpha\|_2^2$$

with $c_{\alpha} \ge 0$ a coefficient of regularization. For example, $b_{\alpha}^s$ can be a polynomial interpolation like in [17] or a spline interpolation where $\alpha$ stands for the coefficients in the polynomial or spline basis. The alternative we choose in this work was to define the surrogate model $b_{\alpha}^s$ as a multilayer perceptron used as a universal approximator like in [3] and in [16]. In that case $\alpha$ stands for the neural network coefficients and they are computed by neural network training. After training, thanks to $b_{\alpha}^s$, we can estimate the value of $b$ at any points $x \in \Omega$ and at any time $t \in [0, T]$; $b_{\alpha}^s(x, t) \approx b(x, t)$.

### Differentiation

The differentiation step gives tools to estimate values of the partial derivatives of $b$ at any point in order to build a catalogue of possible candidate terms for $F_{\beta}$. Obviously, the computational method for the differentiation is related to the nature of the surrogate model. For example, in [17] the extrapolation and hence the partial derivative estimates are performed through a polynomial interpolation technique. In our case, the surrogate model is based on neural networks and the partial derivatives are obtained by back-propagation or automatic differentiation like in [3] and in [16]. Let us denote by $\hat{F}_{\beta}(b_{\alpha}^s(x, t))$ an approximation of $F_{\beta}(b(x, t))$ the right-hand side of Eq. (1) computed at point $(x, t)$ thanks to the surrogate model $b_{\alpha}^s$.

### Identification

All the ingredients are almost there to perform the identification step which consists in solving the following optimization problem. Let us consider $(x_j, t_j)_{j=1,q}$, a set of $q > 0$ points in $\bar{\Omega} \times [0, T]$, not necessarily the same points as the $p$ data points. Eq. (1) is satisfied at any of those points, and we can find the parameters $\beta$ that best fit the data by minimizing the following loss function:
$$J_2(\beta) = \frac{1}{2} \sum_{j=1}^{q} \left| \frac{\partial b_{\alpha}^{s}(x_j, t_j)}{\partial t} - \hat{F}_{\beta}(b_{\alpha}^{s}(x_j, t_j)) \right|^2 + c_{\beta} |\beta|_{1}^{2}$$

where $\frac{\partial b_{\alpha}^{s}(x_j, t_j)}{\partial t}$ is an approximation of the time derivative of b at $(x_j, t_j)$ computed from the surrogate model $b_{\alpha}^{s}$ and $c_{\beta} > 0$ is a real parameter. The second term with 1-norm helps to preserve the sparsity of $\beta$ when there are many candidate terms in the searched PDE. For our work, we chose to implement this optimization problem thanks to a perceptron having the candidate terms as entries and the time partial derivative as output. Forcing the bias to zero, the values of the perceptron coefficients after training corresponds to the values of the parameters $\beta$, ensuring that the resulting PDE is human readable. Fig. 5 summarizes the four steps of the workflow applied to tree biomass, focusing on the data that travel from one step to the next.

Many expert choices must be made to perform the calculation described above: the neural network architecture and hyper-parameters, the design of the loss function, the algorithm to minimize the loss function, the method to compute the partial derivatives for examples. These choices are not easy to make: results obtained with different methods might not be the same. And it still raises theoretical questions and implementation issues [9].

## Method validation

There is no direct way to validate the outputs of the workflow since the governing equations for plant biomass dynamics are not explicitly known. Instead, we monitored the results' validity at each step of the workflow through convergence analysis and consistency with biomass dynamics. In the following, our purpose is 1) to illustrate that the outputs of the workflow are in accordance with the qualitative behavior of the biomass dynamics and 2) to show the workflow can provide quantitative metrics to characterize the dynamics of plants' biomass.

### Experimental settings

We implemented the CEDI workflow in a Python3 script using Tensorflow 2.17 and Keras 3.5. We ran the workflow on CPUs of the MESO@LR-Platform at the University of Montpellier. We used datasets showing the dynamics of three trees of contrasted forms: Massart, Prevost and Rauh architectural types [12]. For each tree, we tested the workflow with four different catalogues or, in other words, four distinct sets of candidate terms for the PDE. We also run the computations with different batch sizes for training. We monitored the loss function during training, we visualized the predicted values of the biomass and compared the obtained PDE models.

### The data

The datasets taken from [1] gather tree biomass values on cubic meshes with step sizes $\delta x_1 = \delta x_2 = \delta x_3 = 20$ cm from zero to 10 years after seedling with time step $\delta t = 0.1$ year. We built these voxelized data from 3D mockups simulated with the AmapSim software [2]. Examples of these mockups are shown on Fig. 2 for the Massart tree, Fig. 3 for the Prevost tree and Fig. 4 for the Rauh tree at t=2, 4, 6, 8 and 10 years after seedling. To ensure the convergence of the training, the dataset has been normalized and scaled according to the following transformations. The normalized biomass has been computed such that $B = \frac{b}{\sigma_b}$ where b is the biomass in gram (g) and $\sigma_b$ is the standard deviation of the biomass values. The time and space coordinates have been scaled as follows: $\tau = \frac{t - \bar{t}}{\sigma_t}$ and $\xi_i = \frac{x_i - \bar{x}_i}{\sigma_i}$, $i = 1, 2, 3$ where $\bar{t}$ and $\bar{x}_i$ are the mean values of time and space coordinates and $\sigma_t$ and $\sigma_i$ are their standard deviations.

### Details of the implementation

The surrogate model of the extrapolation step is a feed-forward neural network made of sequential layers starting with an input layer of four entries for $x_1, x_2, x_3$ and $t$ followed by a 16-node dense layer and then three 32-node dense layers and finally the last layer with one output standing for the biomass. All the internal layers have hyperbolic tangent activation functions and L2-kernel regularization with coefficient set to $10^{-6}$. Preliminary tests showed that higher values (between $10^{-5}$ and $10^{-3}$) tend to erase the spatial and temporal variability of the predicted biomass field, and values higher than $10^{-2}$ fail to converge when training the surrogate model. The choice of hyperbolic tangent is motivated by regularity assumptions required for the differentiation step. The training was

<figure><img src="image_8.png" /><figcaption>Fig. 5. Data flowchart summarizing the CEDI workflow applied to tree biomass dynamics.</figcaption></figure>
performed with Adam optimizer with learning rate set to $10^{-4}$ and Keras default values for other parameters. A learning rate higher than $10^{-3}$ gave unstable training in preliminary tests, while lower than $10^{-5}$ led to slow convergence. We trained the network with batch size equal to 2, 16 and 64 and observed the results after 100 epochs. In the derivation step, we computed estimates of first and second order derivatives of *B* with the Tensorflow function GradientTape using automatic differentiation technique. We computed the derivatives of *B* on the same space and time points as those of the dataset. And for the identification step, we build 4 different catalogues, numbered from I to IV, of candidate terms given in Table 2. These terms are chosen in accordance with available expertise on plant biomass dynamics, see for example [5] and [4]. We assume that a tree may be measured through its global biomass *B* that varies along time and space due to tree growth in height and diameter (growth speed may change over time and direction). These variations depend on the tree species. Translated to mathematical terms, we constructed 4 operator catalogs based on the biomass *B* and its partial derivatives expressed either with Cartesian coordinates or cylindrical coordinates. The latter may offer a good fit since the trees represented in the dataset are roughly axially symmetrical around $\vec{e}_3$ axis. The term *B* is associated with exponential growth, the first order derivatives are associated with transport phenomena and the sum of the second order derivatives is associated with isotropic growth speed of the tree crown. Each catalogue lists a different combination of possible candidate terms for the PDE to be discovered. Catalog I fits with exponential biomass growth, vertical growth velocity and isotropic growth speed of the crown. Catalog II fits with exponential biomass growth, vertical/lateral growth velocity and isotropic growth speed of the crown. Catalog III fits with exponential biomass growth, vertical and radial growth velocity and crown isotropic growth speed in cylindrical coordinates. Catalog IV is equivalent to catalog II with cylindrical coordinates.

Each catalogue gives rise to a different perceptron architecture: the inputs of the tested perceptrons are the terms listed in the catalogue. The output is always the time partial derivative. Again, these data are normalized to ensure a better convergence of the perceptron training. Let us denote *Y* the vector containing the values of $\partial B / \partial t$ (that is the perceptron's output) and $X_i$ the vector containing the $i^{th}$ candidate (that is the $i^{th}$ perceptron's input). The normalization are performed according to the following transformation: $Y^* = \frac{Y}{\sigma_Y}$ and $X_i^* = \frac{X_i}{\sigma_X}$ where $\sigma_Y$ and $\sigma_X$ are the standard deviation of *Y* and $X_i$ respectively. The perceptron was trained with Adam optimizer with learning rate set to $10^{-6}$ and batch size set to 2, 16 and 64. The value of the learning rate could have been increased, since values lower than or equal to $10^{-6}$ lead to slow convergence. A L1-regularization term was set to $10^{-2}$ for sparsity of the weights and the bias was forced to zero. Setting the L1-regularization term below $10^{-3}$ could have been a better choice, since this hyperparameter tends to nullify the values of the searched weights for values above $10^{-2}$.

## Experimental results

### Convergence analysis and stability

Fig. 6 shows the convergence histories obtained when training the surrogate model (step Extrapolation). Convergence is almost reached after 100 epochs. For the Massart tree the smallest loss function value is obtained with a batch size of 16, while the biggest values are obtained with the Prevost tree for a batch size of 64. For the Rauh tree the 3 different values of the batch size seem to converge toward the same values of loss function. For the identification step performed with the catalogues of Table 2, we have seen on Fig. 11 that the batch size has a little impact on the value of the loss function. The batch size seems to impact only the speed of convergence. The learning rate which was set to $10^{-6}$ gave slow convergence: values between $10^{-5}$ and $10^{-3}$ could have led to faster convergence. Catalogues I and II, which are built with derivatives written in the Cartesian coordinate system, gave rather high values of the loss function, reaching plateau at values between 0.5 and 1.5. With the cylindrical coordinate system (catalogues III and IV), the plateaus are slightly lower than in the Cartesian case. However, the batch size impacts the values of the coefficients (or weights) computed during the identification step as shown on Table 3. For example with catalogue number II, from one value of the batch size to the other, the values of $v_1$ and $v_2$ can be totally different. This can be due to several reasons related to the data themselves or the algorithm. The space and time steps of the data might not be fine enough to generate a better solution space and the batch size and batch contents obtained after training data shuffling might not capture some data homogeneity. The data structure, or the loss function, or the starting guess of the optimizer, or the conjunction of all these reasons result in a solution space with many local minima. The same observation can be done for $v_\theta$ in catalogue IV, while catalogue III leads to the convergence of $v_r$ and seems to be a better choice in terms of solution space. A comprehensive study examining the influence of batch size and optimization settings on PDE identification was deferred to future work, though such analysis would greatly enhance the reproducibility and robustness of our

<table><caption>**Table 2**
List of the candidate terms for the 4 catalogues tested.</caption><thead><tr><th></th><th>Candidate terms for the PDE</th><th>Labels of the associated weights in Table 3</th></tr></thead><tbody><tr><td>Catalogue I</td><td>B: <sup>∂B</sup>/<sub>∂x<sub>3</sub></sub>; <sup>∂<sup>2</sup>B</sup>/<sub>∂x<sub>1</sub><sup>2</sup></sub> + <sup>∂<sup>2</sup>B</sup>/<sub>∂x<sub>3</sub><sup>2</sup></sub> + <sup>∂<sup>2</sup>B</sup>/<sub>∂x<sub>3</sub><sup>2</sup></sub></td><td>ρ; v<sub>3</sub>; d</td></tr><tr><td>Catalogue II</td><td>B: <sup>∂B</sup>/<sub>∂x<sub>1</sub></sub>; <sup>∂B</sup>/<sub>∂x<sub>2</sub></sub>; <sup>∂<sup>2</sup>B</sup>/<sub>∂x<sub>3</sub><sup>2</sup></sub> + <sup>∂<sup>2</sup>B</sup>/<sub>∂x<sub>2</sub><sup>2</sup></sub> + <sup>∂<sup>2</sup>B</sup>/<sub>∂x<sub>3</sub><sup>2</sup></sub></td><td>ρ; v<sub>1</sub>; v<sub>2</sub>; v<sub>3</sub>; d</td></tr><tr><td>Catalogue III</td><td>B: <sup>∂B</sup>/<sub>∂r</sub>; <sup>∂B</sup>/<sub>∂x<sub>3</sub></sub>; <sup>∂<sup>2</sup>B</sup>/<sub>∂x<sub>1</sub><sup>2</sup></sub> + <sup>∂<sup>2</sup>B</sup>/<sub>∂x<sub>2</sub><sup>2</sup></sub> + <sup>∂<sup>2</sup>B</sup>/<sub>∂x<sub>3</sub><sup>2</sup></sub></td><td>ρ; v<sub>r</sub>; v<sub>3</sub>; d</td></tr><tr><td>Catalogue IV</td><td>B: <sup>∂B</sup>/<sub>∂r</sub>; <sup>1</sup>/<sub>r</sub> <sup>∂B</sup>/<sub>∂x<sub>3</sub></sub>; <sup>∂<sup>2</sup>B</sup>/<sub>∂x<sub>1</sub><sup>2</sup></sub> + <sup>∂<sup>2</sup>B</sup>/<sub>∂x<sub>2</sub><sup>2</sup></sub> + <sup>∂<sup>2</sup>B</sup>/<sub>∂x<sub>3</sub><sup>2</sup></sub> + <sup>∂<sup>2</sup>B</sup>/<sub>∂x<sub>3</sub><sup>2</sup></sub></td><td>ρ; v<sub>r</sub>; v<sub>θ</sub>; v<sub>3</sub>; d</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>Fig. 6. Convergence history of the extrapolation step for the data of Massart, Rauh and Prevost trees and batch size (BS) 2, 16 and 64.</figcaption></figure>

<table><caption>**Table 3**

Results of the Identification step for catalogues (Cat) I to IV.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Batch size 2</th><th colspan="3">Batch size 16</th><th colspan="3">Batch size 64</th></tr><tr><th>Massart</th><th>Prevost</th><th>Rauh</th><th>Massart</th><th>Prévost</th><th>Rauh</th><th>Massart</th><th>Prévost</th><th>Rauh</th></tr></thead><tbody><tr><td>Cat I</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ρ</td><td>0.171513</td><td>0.125952</td><td>0.108300</td><td>0.144744</td><td>0.356281</td><td>0.123811</td><td>0.182230</td><td>0.557582</td><td>0.180562</td></tr><tr><td>ν<sub>3</sub></td><td>1.636139</td><td>7.886762</td><td>12.502697</td><td>2.081908</td><td>14.352541</td><td>16.067635</td><td>1.057330</td><td>23.830688</td><td>15.047523</td></tr><tr><td>d</td><td>0.000007</td><td>0.000001</td><td>0.000002</td><td>0.000036</td><td>0.012792</td><td>0.000114</td><td>0.003363</td><td>0.038559</td><td>0.005950</td></tr><tr><td>Cat II</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ρ</td><td>0.172407</td><td>0.124174</td><td>0.105447</td><td>0.144677</td><td>0.435770</td><td>0.123829</td><td>0.183556</td><td>0.659764</td><td>0.199495</td></tr><tr><td>ν<sub>1</sub></td><td>-0.007371</td><td>-0.030532</td><td>0.015476</td><td>0.003161</td><td>0.097990</td><td>-0.026351</td><td>-0.009714</td><td>0.240100</td><td>0.036013</td></tr><tr><td>ν<sub>2</sub></td><td>1.159408</td><td>-3.257215</td><td>-2.676747</td><td>0.816525</td><td>6.716690</td><td>0.142282</td><td>1.438752</td><td>13.936645</td><td>3.062459</td></tr><tr><td>ν<sub>3</sub></td><td>1.665916</td><td>7.781262</td><td>11.964144</td><td>2.046030</td><td>19.217773</td><td>15.933753</td><td>1.228995</td><td>27.171227</td><td>16.209906</td></tr><tr><td>d</td><td>0.000005</td><td>0.000010</td><td>0.000005</td><td>0.000038</td><td>0.015663</td><td>0.000099</td><td>0.003496</td><td>0.043053</td><td>0.007138</td></tr><tr><td>Cat III</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ρ</td><td>0.023123</td><td>0.043417</td><td>-0.000052</td><td>0.000280</td><td>0.132076</td><td>0.001193</td><td>0.013196</td><td>0.404871</td><td>0.038897</td></tr><tr><td>ν<sub>r</sub></td><td>0.255366</td><td>0.295080</td><td>0.161790</td><td>0.255238</td><td>0.316948</td><td>0.162605</td><td>0.276321</td><td>0.325743</td><td>0.194013</td></tr><tr><td>ν<sub>3</sub></td><td>-2.018375</td><td>7.817176</td><td>10.520422</td><td>-0.921041</td><td>9.668255</td><td>13.520478</td><td>-2.552996</td><td>17.761522</td><td>12.505530</td></tr><tr><td>d</td><td>-0.000000</td><td>-0.000000</td><td>0.000000</td><td>-0.000073</td><td>0.002752</td><td>-0.000071</td><td>0.000342</td><td>0.017509</td><td>0.003785</td></tr><tr><td>Cat IV</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ρ</td><td>0.021942</td><td>0.042844</td><td>-0.002183</td><td>0.000322</td><td>0.144933</td><td>0.001171</td><td>0.014234</td><td>0.440772</td><td>0.062709</td></tr><tr><td>ν<sub>r</sub></td><td>0.256145</td><td>0.292236</td><td>0.162157</td><td>0.255994</td><td>0.312831</td><td>0.168146</td><td>0.275941</td><td>0.305816</td><td>0.180212</td></tr><tr><td>ν<sub>θ</sub></td><td>-0.470653</td><td>0.532653</td><td>1.010458</td><td>-0.209571</td><td>2.873393</td><td>-0.866930</td><td>-0.135464</td><td>6.037452</td><td>0.577699</td></tr><tr><td>ν<sub>3</sub></td><td>-2.080043</td><td>7.960285</td><td>10.107252</td><td>-0.971721</td><td>10.907381</td><td>13.683924</td><td>-2.532960</td><td>21.319525</td><td>12.272537</td></tr><tr><td>d</td><td>-0.000003</td><td>0.000002</td><td>-0.000000</td><td>-0.000075</td><td>0.003444</td><td>-0.000065</td><td>0.000390</td><td>0.018244</td><td>0.004817</td></tr></tbody></table>

study.

## Validity of the results with respect to biomass dynamics

The values of the tree biomass predicted by the surrogate models (step Extrapolation) inside the plane $(O, \vec{e}_2, \vec{e}_3)$ shown on Fig. 7, Fig. 8 and Fig. 9 are visually in accordance with the data. Besides, the mean absolute errors between the predicted scaled biomass and the scaled data were equal to 0.13 for Rauh and Massart, while they were equal to 0.26 for Prevost with the three different batch sizes. Some of the partial derivatives computed by automatic differentiation are shown on Fig. 10 for the Massart tree. Again, the plots show the values inside the $(O, \vec{e}_2, \vec{e}_3)$-plane. Visually, the first order derivatives are consistent with the biomass dynamics. The second order partial derivatives seem to be blurred by numerical errors which may affect the rest of the workflow. The computed parameters can be interpreted with respect to the biomass dynamics. For all the trees and all the catalogues tested, the growth rates $\rho$ associated with exponential growth of the biomass is rather small, ranging between almost zero and 0.6 per year. Exponential growth might not be a good choice for these datasets, instead we could have considered logistic growth for example. Similarly, the diffusion $d$ which is associated with the second order derivatives is really small. This might be due to the fact that either diffusion occurs with really small amplitude or the second order derivatives are poorly approximated (see Fig. 10). Additional runs should be performed with
<figure><img src="image_10.png" /><figcaption>Fig. 7. Biomass of the Massart tree plotted on the ($O, \vec{e}_2, \vec{e}_3$)-plane: data (top) and predicted values (bottom) by the extrapolation with the surrogate model (trained with a batch size of 2).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Biomass of the Prevost tree plotted on the ($O, \vec{e}_2, \vec{e}_3$)-plane: data (top) and predicted values (bottom) by the extrapolation with the surrogate model (trained with a batch size of 2).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 9. Biomass of the Rauh tree plotted on the ($O, \vec{e}_2, \vec{e}_3$)-plane: data (top) and predicted values (bottom) by the extrapolation with the surrogate model (trained with a batch size of 2).</figcaption></figure>

different hyperparameter values, in particular with a smaller L1-regularization term, to verify whether the diffusion still cancels. The transport of the biomass along the vertical axis ($v_3$) also called vertical velocity is well captured for the Prevost and Rauh trees. The values of $v_3$ are higher for Rauh than for Prevost, which is in accordance with the fact that the Rauh tree reached higher altitude than Prevost after 10 years of growth: 5 meters against less than 2 meters. With these values, we would have expect a vertical velocity ($v_3$) of around 50 centimeters per year for Rauh and 20 centimeters per year for Prevost in the case of continuous growth in time, while the computations gave smaller values of $v_3$ for all the trees. These differences might be due to the fact that the transport term captures variations in leaf biomass and leaf turnover along the life of the trees. Catalogues I and III gave better results (in particular for the
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Massart tree. Plots on the ($O, \vec{e}_2, \vec{e}_3$)-plane of isovalues of the biomass partial derivatives associated with catalogue II. The first line shows isovalues of $\frac{\partial B}{\partial x_1}$ at t=2, 4, 6, 8 and 10 years (from left to right). The following lines show isovalues of $\frac{\partial B}{\partial x_2}$, $\frac{\partial B}{\partial x_3}$ and $\frac{\partial^2 B}{\partial x_1^2} + \frac{\partial^2 B}{\partial x_2^2} + \frac{\partial^2 B}{\partial x_3^2}$ respectively at t=2, 4, 6, 8 and 10 years (from left to right).</figcaption></figure>

transport parameters) than catalogues II and IV respectively in terms of stability. This can be interpreted at the biomass dynamics level by noticing that lateral transport of the biomass occurs roughly, first along the vertical axis and second, along all radial directions with approximately the same amplitude.

## Limitations

Each of the four steps of the CEDI workflow raises modeling questions as well as theoretical and technical issues that we did not address in the paper but that must be considered for future works. Since every step depends on the previous one, the final result may be influenced by some error propagation and modeling choices. Given the CEDI workflow, these issues may probably strongly depend on the PDE that has to be discovered and the kind of subject under study. The dataset processed during the first step (C) must be a relevant representation of the dynamics of interest. A trade-off between representativity, cost of acquisition/curation and accuracy of the data must be handled. In our work, even though the data space-time sampling gives a good representation of the dynamics of the biomass (since they are generated thanks to a dedicated software) we encountered some stability issues that could be worse in cases of data collected in the field. Acquiring data with such small space and time sampling steps can be a tedious and costly task, unless one uses remote sensing tools for example. The challenge with the design of the surrogate model during the second step (E) is then to compute missing data that is in line with the dynamics under study. In addition, the setting of the surrogate model raises issues about the choice and the implementation of an extrapolation method, which is expected to give a faithful representation of the dynamics at any space-time point. Computing estimates of partial derivatives (step D) requires that the surrogate model is differentiable, which has not been proved here (and which is not necessarily always the case). With the neural network approach, we rely on the automatic differentiation tools implemented in deep learning libraries, but we did not estimate the numerical errors made. The identification step (I) also raises theoretical and implementation issues. This kind of problem is often ill-posed and convergence and stability may be analyzed carefully. Implementation issues particularly deal with the management of boundary conditions, the design of a neural network architecture and the setting of a successful optimizer. As each step depends on the previous one, an accumulation of errors can occur, which can lead to a
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Convergence histories of the identification step with catalogues I to IV for Massart, Prevost and Rauh trees with batch sizes (BS) 2, 16 and 64.</figcaption></figure>

lack of stability or bad convergence rate in the identification step. Studies need to be conducted to evaluate the error propagation behavior inside the CEDI workflow. For each specific case, this may be achieved by the collaboration of domain experts who provide multiple empirical datasets with analytical point of view and numerical experts who perform successive trials and tests on various settings along the four CEDI steps.

Through a simple numerical implementation, we showed that the CEDI workflow is able to capture the qualitative dynamical behaviour of trees' biomass. However, our attempt for the modeling of tree dynamics is far from being perfect in view of the architecture but also the physiology of the plants. The work presented in the paper has only the merit of illustrating the workflow in a concrete case, as proof of concept and to shed light on the many difficulties encountered and improvements to be done. The dynamics of plants are complex in many and various aspects. Plants are living and growing organisms that can be studied at different times and spatial scales. In addition to all the intraspecific and interspecific variability linked to the plants themselves, the environment has a strong influence on plant growth and dynamics through soil conditions, water availability, air and soil temperatures, human practices, pests and diseases for example. Plant growth and dynamics are driven by many physiological, mechanistic, hydraulic and agronomic processes. Few modeling approaches link all these processes in a single model. Each process is associated with its specific type of data which has been massively collected for a very long time, making modeling even more complex, and integrating all that knowledge is a complex task. The CEDI workflow offers a generic tool that might open new roads for the characterization of the plant biomass dynamics provided meaningful datasets are available. We also are convinced that this genericity has to be carefully controlled and set up through a collaboration between experts on the current topic and experts on numerical studies.

To go further, future work concerns the application of the CEDI workflow to real-world field data, since the design of the workflow was originally motivated by the characterization of biomass dynamics in cocoa-based agroforestry systems. Remote sensing technologies such as airborne or terrestrial Lidar provide tools to monitor ecosystems at different scales. Lidar data of real trees consists of 3D point clouds (see for example the dataset described in [15]) whose use in the CEDI workflow could raise many challenges, particularly in the data Collection step. First, a Lidar scan provides static information at a single date while the life of a tree may last several decades. The monitoring of tree biomass dynamics along its life requires many scans performed at well chosen dates according to the phenology of the tree, which would lead to long and expensive data collection campaigns. An alternative based on a strong knowledge of the phenology would consist in reconstructing chrono-sequences of Lidar data from scans of different trees at different stages of their life, provided they all have grown in similar conditions. Another alternative would consist in morpho-architectural studies of scans of a single tree to reconstruct its past. Second, even though Lidar scans provide tree descriptions with high space
accuracy, they suffer from artefacts such as occlusions caused by hidden vegetation or blurring due to winds. To limit occlusions, several scans have to be done at different positions around the same single tree, increasing the operating time in the field and the memory size of the data. To limit the blurring effect, the tree branches and leaves must not move between two scanning positions, requiring good weather conditions without wind nor rain. Third, as mentioned above Lidar data consists of 3D point clouds characterized by the list of point coordinates that have been hit by the laser ray of the Lidar scanner device. The development of algorithms dedicated to the analysis of such data is a research question in itself where the challenge is to assess traits such as woody biomass and foliage distribution from point cloud coordinates. Fourth, to tackle realistic applications that have a real interest for plant scientists, additional traits on top of the tree biomass must be taken into account since issues in the management of agro-ecosystems deal with yield, light availability, pest and disease control, climate change mitigation among others. These issues require applying the CEDI workflow not only at the tree scale but at the plot scale, which involves manipulating bigger and more complex datasets and leading to challenges linked to scaling up, potentially noisy and incomplete data. Finally, many back-and-forth loops between plant and digital experts would be needed to build meaningful dataset in relation to the scientific question and capture the dynamics with sufficient representativeness to enable stable and accurate numerical processing of steps E, D and I of the workflow.

## Conclusion

The main contribution of our work is to introduce a generic workflow for the complex task of discovering PDE models from data. Our workflow, called CEDI, is divided into four main steps: the data Collection, Extrapolation, Derivation and identification steps. The four steps encompass the task of gathering a meaningful dataset, the setting up of a surrogate model in case of missing data, the computation of partial derivatives and finally the extraction of the PDE that best matches with the dataset. The workflow CEDI offers a conceptual framework that unifies methodologies for extracting PDE from data, and has the advantage of being generic and flexible. The workflow is particularly useful for applications where theoretical models are incomplete or unavailable. As a preliminary candidate, we applied the workflow to the modeling of tree biomass dynamics using realistic data and a neural network approach. While this basic application gave promising results, future work should focus on improving the neural network architectures and a comprehensive study of the impact of the hyperparameters on the convergence and accuracy. The numerical results presented in this paper have the merit of illustrating the conceptual framework of the CEDI workflow and is a methodological milestone toward applications to real world data, to be performed in conjunction with plant and digital experts, in the perspective to address questions related to the understanding of plants' dynamics. The workflow paves the way for cognitive tools to generate knowledge through data-driven PDE modeling within the paradigm of theory-driven data science, while enabling researchers from diverse fields to pool their efforts and structure their work in favor of the understanding of dynamical phenomenon in life sciences.

## Ethics statements

Our work did not involve human subjects, animal experiments or data collected from social media platforms.

## CRediT author statement

**Emilie Peynaud:** Conceptualization, Methodology, Software, Validation, Formal analysis, Writing - Original Draft, Writing - Review & Editing, Visualization, Supervision, Project administration, Funding acquisition. **Paulin Melatagia:** Conceptualization, Formal analysis, Methodology, Validation, Writing - Original Draft, Writing - Review & Editing. **Serge Stinckwich:** Conceptualization, Methodology, Writing - Original Draft, Writing - Review & Editing. **Jean-François Barczi:** Data production, Writing - Review & Editing, Resources, Data Curation, Visualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

This work was supported by Cirad (French Agricultural Research Center for International Development) [grant AI CRESI 2019] and it has been realized with the support of MESO@LR-Platform at the University of Montpellier, France.

## Data availability

Links to data and source code are given in the specifications table.

## References

[1] J.-F. Barczi, Plant Growth Simulation Data, CIRAD Dataverse V1, (2025), https://doi.org/10.18167/DVN1/VK0J16.
[2] J.-F. Barczi, H. Rey, Y. Caraglio, P. de Reffye, D. Barthélémy, Q. Dong, T. Fourcaud, AmapSim: A structural whole-plant simulator based on botanical knowledge and designed to host external functional models, Ann. Bot. 101 (8) (2007) 1125–1138, https://doi.org/10.1093/aob/mcm194.

[3] J. Berg, K. Nystrom, Data-driven discovery of PDEs in complex datasets, J. Comput. Phys. 384 (2019) 239–252, https://doi.org/10.1016/j.jcp.2019.01.036.

[4] R. Beyer, O. Etard, P.-H. Cournède, P. Laurent-Gengoux, Modeling spatial competition for light in plant populations with the porous medium equation, J. Math. Biol. 70 (2015) 533–547, https://doi.org/10.1007/s00285-014-0763-1.

[5] A. Bonneu, Y. Dumont, H. Rey, C. Jourdan, T. Fourcaud, A minimal continuous model for simulating growth and development of plant root systems, Plant Soil. 354 (2012) 211–227, https://doi.org/10.1007/s11104-011-1057-7.

[6] G.J. Both, S. Choudhury, P. Sens, R. Kusters, DeepMoD: deep learning for model discovery in noisy data, J. Comput. Phys. 428 (2021) 109985, https://doi.org/10.1016/j.jcp.2020.109985.

[7] S.L. Brunton, J.N. Kutz, Promising directions of machine learning for partial differential equations, Nat. Comput. Sci. 4 (2024) 483–494, https://doi.org/10.1038/s43588-024-00643-2.

[8] D. Barthélémy, Y. Caraglio, Plant architecture: a dynamic multilevel and comprehensive approach to plant form, structure and ontogeny, Ann. Bot. 99 (2007) 375–407, https://doi.org/10.1093/aob/mcl260.

[9] S. Cuomo, V. Di Cola, F. Giampaolo, G. Rozza, M. Raissi, P. Piccialli, Scientific machine learning through physics-Informed neural networks: where we are and what's next, J. Sci. Comput. 92 (3) (2022) 88, https://doi.org/10.1007/s10915-022-01939-z.

[10] E. Grenier, M. Hoffmann, T. Lelièvre, V. Louvet, C. Prieur, N. Rachdi, P. Vigneaux, Statistical inference for partial differential equations, ESAIM: Proc. Surv. 45 (2014) 178–188, https://doi.org/10.1051/proc/201445018.

[11] G. Frasso, J. Jaeger, P. Lambert, Parameter estimation and inference in dynamic systems described by linear partial differential equations, AStA Adv. Stat. Anal. 100 (3) (2016) 259–287, https://doi.org/10.1007/s10182-015-0257-5.

[12] F. Hallé, R.A. Oldeman, P.B. Tomlinson, Tropical Trees and Forests: An Architectural Analysis, Springer, Berlin, Heidelberg, 1978, https://doi.org/10.1007/978-3-642-81190-6.

[13] A. Karpatne, G. Alturi, J.H. Faghmous, A. Banerjee, A. Ganguly, S. Shekhar, N. Samatova, V. Kumar, Theory-guided data science: A new paradigm for scientific discovery from data, IEEe Trans. Knowl. Data Eng. 29 (10) (2017) 2318–2331, https://doi.org/10.1109/TKDE.2017.2720168.

[14] D. Kim, J. Lee, A review of physics informed neural networks for multiscale analysis and inverse problems, Multiscale Sci. Eng. 6 (2024) 1–11, https://doi.org/10.1007/s42493-024-00106-w.

[15] E. Peynaud, S. Momo Takoudjou, Terrestrial LiDAR point cloud dataset of cocoa trees grown in agroforestry systems in Cameroon, Data Br. 53 (2024) 110108, https://doi.org/10.1016/j.dib.2024.110108, 9 p.

[16] M. Raissi, P. Perdikaris, G. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, J. Comput. Phys. 378 (2019) 686–707, https://doi.org/10.1016/j.jcp.2018.10.045.

[17] S. Rudy, S. Brunton, J. Proctor, J. Kutz, Data-driven discovery of partial differential equations, Sci. Adv. 3 (4) (2017) e1602614, https://doi.org/10.1126/sciadv.1602614.

[18] H. Schaeffer, Learning partial differential equations via data discovery and sparse optimization, Proc. R. Soc. A: Math. Phys. Eng. Sci. 473 (2197) (2017) 20160446, https://doi.org/10.1098/rspa.2016.0446.

[19] R. Stephany, E. Christopher, PDE-READ: Human-readable partial differential equation discovery using deep learning, Neural Netw. 154 (2022) 360–382, https://doi.org/10.1016/j.neunet.2022.07.008.