# Unraveling biochemical spatial patterns: machine learning approaches to the
inverse problem of Turing patterns

Antonio Matas-Gil1, Robert G. Endres1†

(1Department of Life Sciences & Centre for Integrative Systems Biology and
Bioinformatics, Imperial College London, London SW7 2BU, United Kingdom
<EMAIL_ADDRESS>
###### Abstract

The diffusion-driven Turing instability is a potential mechanism for spatial
pattern formation in numerous biological and chemical systems. However,
engineering these patterns and demonstrating that they are produced by this
mechanism is challenging. To address this, we aim to solve the inverse problem
in artificial and experimental Turing patterns. This task is challenging since
high levels of noise corrupt the patterns and slight changes in initial
conditions can lead to different patterns. We used both least squares to
explore the problem and physics-informed neural networks to build a noise-
robust method. We elucidate the functionality of our network in scenarios
mimicking biological noise levels and showcase its application through a
prototype involving an experimentally obtained chemical pattern. The findings
reveal the significant promise of machine learning in steering the creation of
synthetic patterns in bioengineering, thereby advancing our grasp of
morphological intricacies within biological systems while acknowledging
existing limitations.

## Introduction

Spatial patterns are prevalent in biological systems, including gene
expression in microbial communities, developing embryos, as well as skin and
fur patterns in adult animals. A leading mechanism for pattern formation is
the diffusion-driven instability in reaction-diffusion models, as proposed by
Alan Turing [1] and extended by others [2, 3, 4, 5, 6]. These models typically
describe interacting and diffusing activator and inhibitor chemical species
through sets of coupled partial-differential equations (PDEs), although
simulation methods are also available. While Turing patterns were originally
observed in chemical systems [7, 8] it has been challenging to conclusively
demonstrate this mechanism in developmental systems such as digit formation
[9], zebra-fish skin pigment patterning [10] and hair spacing in mice [11]
(for recent reviews see [12, 13]). Additionally, building Turing patterns from
bottom-up with synthetic circuits has proven difficult [14, 15]. Identified
issues are that the Turing mechanism is exceedingly simple compared to
biological regulatory pathways, highly sensitive to changes in model
parameters [16], and experimental control over parameter tuning is limited.
Even if patterns can be produced experimentally, the challenge remains of
linking them to the corresponding parameters of a candidate model to support
it as an actual Turing pattern rather than an experimental artifact. Further
complications arise from the relatively limited amount of data available in
developmental systems, which is often corrupted by measurement and imaging
noise, as well as the strong pattern variability observed in microbial
systems. Addressing these issues would greatly benefit from the ability to
robustly estimate parameters from given patterns and candidate models.

To address the issues that arise in an experimental setting, the inverse
problem can be initially approached using artificial data from numerical
simulations instead of actual experimental data. This involves generating
random initial conditions and evolving the reaction-diffusion model in space
and time, allowing for a systematic study of the effects of data size and
noise. However, solving the inverse problem remains a challenging task even
with this simplification, as the same model parameters can produce different
patterns, resulting in a many-to-one inversion problem [17]. This occurs
because patterns are highly sensitive to initial conditions, and even slight
alterations can cause noticeable changes in the final pattern obtained. These
changes do not affect the type of pattern obtained, such as spots or stripes,
but do alter the location and shape of the pattern elements. As a result,
previous research has ruled out direct minimization of the mean-squared
displacement between data and model output for parameter fitting, considering
it an ill-posed problem [18]. Instead, focus has been on more advanced
approaches, including Bayesian inference and other statistical tools [18, 19]
and different machine-learning techniques such as support vector machines,
Kernel-based methods, and neural networks [20], as well as optimal control
theory [21]. However, despite some success, these approaches suffer from
requiring thousands of training images [19], sensitivity to noise in the
patterns [20], ad hoc cost or loss functions for quantifying the quality of
fit [18], or the requirement to fix some parameters and have knowledge of
initial conditions and the pattern evolution [21].

As the inverse Turing problem remains largely unresolved, there is a need to
develop new robust approaches, particularly when dealing with small and noisy
data. Although conceptually straightforward, the least squares (LS) approach
has not been extensively explored in this context. This method requires the
definition of a parameter-dependent loss function, which can be the residuals
of the model equations [22]. In contrast, deep learning methods, such as
neural networks, are usually more reliable when it comes to noise, even though
they are computationally expensive to train. A key property that makes neural
networks a promising tool for approaching this problem is the universal
approximation theorem, which states that given enough parameters, neural
networks are capable of approximating any continuous function, including
spatial patterns [23, 24]. Additionally, physics-informed neural networks
(PINNs) incorporate physical constraints, such as the model equations that
should hold for the data, thereby helping to regularize the training [25] and
making PINNs more intuitive than black-box neural networks. An added benefit
of PINNs is that they do not require large amounts of data; in one case, a
single simulation output was sufficient to solve the inverse problem [26].
Model parameters are regarded as parameters of the neural network, resulting
in the computational cost of the whole method being of the same order as the
standalone function approximation. This makes PINNs superb candidates for
solving the inverse problem in Turing patterns.

Here, we explore basic LS applied to the PDE residuals and advanced PINNs to
address the inverse problem - given Turing patterns and candidate models, we
aim to recover the model parameters. Using our first approach, we find that LS
is computationally inexpensive, but that it requires mostly exact (albeit
little) data and hence does not allow significant noise or using a similar
pattern produced from a different model. Applying this methodology to small
regions of patterns still allows us to recover the parameters, and we identify
the minimum number of necessary pixels. Our second approach uses physics-
informed neural networks with a radial basis function (RBF) architecture to
approximate the patterns, referred to as RBF physics-informed neural networks
(RBF-PINNs). This method remedies many of the issues from the first method at
a two-order of magnitude higher computational costs, allowing us to obtain
accurate results up to a 10-20% relative noise given a single snapshot of a
Turing pattern, even in an experimental setting with chemical Turing patterns.
Our least-squares results show a new perspective on the Turing robustness
problem, that a pattern can act like a barcode to a specific parameter
combination, given the key of the correct model. Our RBF-PINNs are a promising
method to solve the inverse problem, and to guide future experiments in
bioengineering in the development of synthetic tissues.

## Results

### Models

Both of our approaches, LS and RBF-PINNs, utilize discrete Turing patterns as
inputs to infer the parameters of a candidate model. In this work, the
candidate model refers to the true model that was numerically solved to
produce the images using the finite-differences algorithm (see Methods). We
explored several models, all following the idea of activator-inhibitor
dynamics shown in Fig. 1A but before delving into specific ones, we present
the general two-component reaction-diffusion model for concentrations $u$ and
$v$, which depend on space and time, as follows:

$\displaystyle\begin{split}u_{t}&=D_{u}\Delta u+f(u,v)\\\ v_{t}&=D_{v}\Delta
v+g(u,v)\end{split}$ (1)

where $u_{t}$ and $v_{t}$ denote partial differentiation with respect to time
and the functions $f(u,v)$ and $g(u,v)$ are non-linear reaction kinetics.
Depending on the specific $f$ and $g$ functions, different models can be
identified. We will focus on three of these models. First, the Schnakenberg
model [2], which has a total of 6 parameters with 4 of them inside the
functions $f$ and $g$ and the other two given by the diffusion coefficients:

$\displaystyle\begin{split}f(u,v)&=c_{1}-c_{2}u+c_{3}u^{2}v,\\\
g(u,v)&=c_{4}-c_{3}u^{2}v\end{split}$ (2)

Second, we have the FitzHugh-Nagumo model [5, 6], which has 5 parameters:

$\displaystyle\begin{split}f(u,v)&=c_{1}(v-c_{2}u),\\\
g(u,v)&=-u+c_{3}v-v^{3}\end{split}$ (3)

Third, the Brusselator model [4], with a total of 4 parameters:

$\displaystyle\begin{split}f(u,v)&=c_{1}-(c_{2}+1)u+u^{2}v,\\\
g(u,v)&=c_{2}u-u^{2}v\end{split}$ (4)

We numerically solve these PDE models on a square domain with zero-flux
boundary conditions, resulting in patterns similar to Fig. 1B. There are
different types of patterns, e.g. dots (Fig. 1B, top left), and labyrinths
(Fig. 1B, top right). The type is mostly dependent on the model that we
choose: FitzHugh-Nagumo produces labyrinths while Schnakenberg produces dots;
but some models, like the Brusselator, can produce several types of patterns
depending on the parameters. Parameters for the Brusselator and the FitzHugh-
Nagumo models can be found in [18], while those for the Schnakenberg model are
in [27]. Different initial conditions will produce different patterns of the
same type, but fixing these will produce the same one. Hence, we can
conceptually think of the parameters of the model and the initial conditions
as the ‘variables’ that produce a given pattern. Since we are mostly
interested in the parameters, we used the same initial conditions for all
patterns. This eases comparison between the patterns produced by different
methods. As a result, there is a more direct relationship between the
parameters and the patterns. The inverse problem consists of inverting this
relationship by recovering the parameters that generate a given pattern, as
outlined in Fig. 1C.

Figure 1: Turing patterns and methods for the inverse problem. (A) Network of
canonical Turing pattern with a short-range activator and a long-range
inhibitor. (B) Turing patterns from different models. The top-left pattern
showing spots is produced with the Schnakenberg model, the top-right pattern
showing labyrinths with FitzHugh-Nagumo and the bottom row, also showing spots
and labyrinths, with the Brusselator and two different parameter sets. (C)
Methodology of the paper. Starting from a Turing patterns, we build a PDE loss
which is minimized with respect to the parameters $\beta$. This is done using
two different methods, LS and RBF-PINNs.

### Non-dimensionalization

Depending on the model, the inverse problem described above may be unsolvable
because it may have many solutions. For example, take the Schnakenberg model,
Eq. 1 with reaction kinetics given by Eq. 2, and assume $u$ and $v$ are
steady-state patterns in time satisfying the PDE (Eq. 1) for a given set of
parameters. As a result, both the left-hand side (LHS) and right-hand side
(RHS) of the PDE are zero. If we multiply all parameters by a constant $k$ and
substitute them back in the equation, the resulting set of parameters is also
a solution of the inverse problem; since we can take $k$ as a common factor,
we again obtain $RHS=0$. Furthermore, if we let $k=0$ we arrive at a trivial
solution where all parameters are zero. To avoid this, one remedy is to fix
some parameters so that there is only a single solution. However, we found
that non-dimensionalizing the equations to decrease the number of parameters
is preferable, since this does not require us to make any assumptions on the
parameters. As an added benefit, this reduces the number of parameters of the
model.

We can write the previous models in non-dimensional form as follows. For the
Schnakenberg model we obtain:

$u_{t}=\Delta u+c_{1}-u+c_{2}u^{2}v\qquad v_{t}=d\Delta
v+c_{3}-c_{2}u^{2}v,\\\ $ (5)

The FitzHugh-Nagumo model can be rewritten as:

$u_{t}=\Delta u+v-c_{1}u\qquad v_{t}=d\Delta v+c_{2}v-c_{3}u-v^{3},\\\ $ (6)

And the non-dimensional form of the Brusselator model is:

$u_{t}=d\Delta u+1-u+c_{1}u^{2}v\qquad v_{t}=\Delta v+c_{2}u-c_{1}u^{2}v$ (7)

where in all equations $d$ is the ratio of the diffusion coefficients and the
rest of parameters are non-dimensional. Note that even though we have not
changed the notation of the parameters, they are different from the ones in
Eqs. 2, 4 and 3. For the definition of the new variables and parameters for
the three models see Supplementary Information section 2.

### Least squares

Our first approach to solve the inverse problem is based on fitting the
parameters to the PDE equations, by fixing the concentrations $u$ and $v$ as
given by the Turing patterns we aim to reproduce. If we assume the pattern to
be a steady state of the PDE in e.g. Eq. 5, then $u_{t}=v_{t}=0$. We can now
consider the RHS of Eq. 5. Because we assume we have access to the patterns
$u,v$ satisfying the PDE, we can fix them and treat the RHS as a function
dependent only on the parameters, which will be zero for the combination of
parameters that we want to find. Since our patterns are discrete in space (we
can think of them as images with pixels), $\mathbf{u}$ and $\mathbf{v}$ are
matrices, so we can write its elements as $u_{ij}$ and $v_{ij}$ for
$i,j=1,2,...,N$, where $N$ is the number of rows and columns in the pattern.
This makes it possible to write and solve the problem using LS, as we can now
formulate the problem as $X\beta=Y$, where $X$ is the design matrix, $\beta$
is the vector of parameters and $Y$ are all terms that are not parameter
dependent. (For an example, see the Methods section.) This LS minimization is
equivalent to minimizing the squared of the Frobenius norm (element-wise
$L_{2}$ norm for matrices) for both equations, which we can write as:

$L(\boldsymbol{\beta})=||\beta_{1}\Delta\mathbf{u}+f(\boldsymbol{\beta})||_{F}^{2}+||\Delta\mathbf{v}+g(\boldsymbol{\beta})||_{F}^{2}$
(8)

Figure 2: Least squares for parameter inference. (A) In LS, noise changes the
loss function that we are minimizing and the optimal value shifts, leading to
some error in the inferred parameters. (B) Resulting patterns obtained from LS
with different noise levels. After corrupting the original pattern with
relative noise of different levels, parameters are obtained using LS and the
model is solved with the inferred parameters to obtain the new patterns shown.
The upper row corresponds to the FitzHugh-Nagumo model, the lower row to the
Schnakenberg model. Enlarged patterns on the left are the original ones. It
can be observed that at 1% noise, the predicted patterns are noticeably
different from the original ones across the models, and that we can obtain
several patterns at the same level of noise. (C) Radially averaged power
spectra (RAPS) obtained for different patterns recovered at different noise
levels (shown in (B)) from the Schnakenberg model. Red corresponds to 1%
noise, green to 0.6%, orange to 0.3% and blue is the original pattern. Inset
shows the same plot but the $y$ axis in log-scale. (D) Scatter plot of RAPS
differences for different relative noise levels with noticeable sudden
increases due to discrete changes in the predicted Turing patterns. (E)
Average relative error of inferred parameters for different relative noise
levels in the Schnakenberg and FitzHugh-Nagumo models, with a sketch of the
parameter space explaining the difference in spread. Green point represents
original set of parameters values and orange points are the different optimal
sets resulting from the minimization. The bigger standard deviation in
relative error occurs when the optimal sets are close to the original so the
scale of the error changes drastically (can be very close or far), while the
smaller standard deviation occurs when the optimal points are further away and
the error is always on a similar scale.

The LS method has previously been regarded in the literature [23], but to our
knowledge it has not been thoroughly investigated. The Laplacian is
approximated using a second-order finite difference method, and the no-flux
boundary conditions are taken to be the same as used to simulate the pattern.
Since we assume the Turing pattern is at steady state, so that $u_{t}=0$, we
refer to this method as a steady-state approximation. The LS method can also
be extended to dynamic patterns, e.g when having two samples of the time
evolution of the system, say $u_{1}$ and $u_{2}$ at $t_{1}$ and $t_{2}$,
respectively. In this case, we can no longer say that $u_{t}=0$, but instead
approximate this partial derivative using finite differences so that we obtain
$u_{t}\approx\frac{u_{2}-u_{1}}{t_{2}-t_{1}}$ (see Discussion).

Figure 3: Effect of the number of pixels in the least squares method. (A)
Pattern from the FitzHugh-Nagumo model with three different cropped regions
given by small squares of different size. It can be seen that a region of
$3\times 3$ pixels is sufficient to recover accurate enough parameters such
that the predicted and original patterns (right and left respectively) are
indistinguishable. (B) Schematic of choosing $N$ randomly selected points
(black) on the Turing pattern. (C) Effect of increasing the number of randomly
selected pixels on the average relative error in the inferred parameters with
(orange) and without (blue) added noise to the original pattern for the
FitzHugh-Nagumo model. We used $N$ in the range $5$ to $2000$ ($50\times
50=2500$ being the maximum possible) and sampled $10$ different sets of pixels
for each $N$, and we measured the relative error for each of the inferred
parameters. There is almost no effect without noise, but with noise there is a
steady reduction in the relative error. Also shown is the slope of the line of
best fit to the data (orange).

To apply the LS method, we first select a model and a set of parameters,
$\beta_{orig}$, that can produce a Turing pattern. Then, using the model and
$\beta_{orig}$ we numerically solve the PDEs to produce a pattern similar to
the ones shown in Fig. 1B. Once we have the patterns, we minimize
$L(\boldsymbol{\beta})$. This method produces very accurate results without
noise in the pattern, with an average relative error in parameters of the
order of $O(10^{-15})$ or lower, which can be considered artifacts given its
closeness to numerical precision. When we add noise to the pattern before
applying LS, the ‘true’ minimum shifts from $\beta_{orig}$, producing an error
which is no longer a numerical artifact (Fig. 2A). The way we incorporate
noise is by adding a matrix of normally distributed random variables with zero
mean and varying standard deviation $\sigma$ to the concentration matrices
$\mathbf{u}$ and $\mathbf{v}$. Since different $\mathbf{u}$ and $\mathbf{v}$
will have different ranges, the standard deviation of the noise must be
relative to the range (maximum minus minimum value in the pattern) of each
concentration matrix. To achieve this we employ what we called ‘relative
noise’. If we let $R_{p}(\mathbf{u})$ be the range of the concentration matrix
$\mathbf{u}$, so that

$R_{p}(\mathbf{u})=\max_{1\leq i,j\leq N}\\{u_{ij}\\}-\min_{1\leq i,j\leq
N}\\{u_{ij}\\}$

we can define a $s\%$ relative noise to correspond to a standard deviation
$\sigma=\frac{s}{100}R_{p}$. Once we add noise to the pattern, we can use
these ‘noisy patterns’ as input for our methods, and by changing the level of
noise we can systematically investigate the difference in performance and
accuracy. This will yield a different parameter set for each noisy pattern,
which can be used to solve the system and obtain a new pattern. The new
patterns obtained and the respective level of noise from which their
parameters were inferred are shown in Fig. 2B. Before delving into our
results, we remark that the figures only show the obtained results for the
Schnakenberg and FitzHugh-Nagumo models. All the corresponding figures for the
Brusselator model can be found in Supplementary Information section 3 and will
be referenced when necessary.

We considered several measures to assess the accuracy of this method. First,
we looked at the average relative error in the parameters, which gives us a
measure of how close the newly obtained parameters are to the original ones. A
drawback of this measure is that, depending on the model, very distinct
parameters can give us the same pattern, or conversely similar parameter
values can yield very different patterns. Hence, a better measure for
parameter accuracy would be one that quantifies how similar the pattern
produced with an inferred parameter vector is to the original pattern. We
found that measures like the mean squared error (MSE) are not very useful,
since we do not expect the patterns to be exactly the same. This is because
our main objective is to be able to recover the parameters of the model, but
the final pattern is not determined solely by the parameters, but also by the
initial conditions; the parameters and model determine the type of Turing
pattern and its wavelength, but the position and shape are determined by the
initial conditions. Hence, we do not put emphasis on obtaining the exact
positions and shapes. In order to focus on the type of pattern and wavelength,
we instead compare the patterns in the frequency domain. Specifically, we use
the Fourier power spectrum and take the radial average to obtain a one-
dimensional profile which should have a main dominant frequency for each of
our patterns. This is called radially averaged power spectrum (RAPS). Example
RAPS curves for different profiles are shown in Fig. 2C. Then, we can compute
this RAPS for each pattern and use the MSE between these two profiles as a
similarity measure. We will refer to this measure as RAPS difference or RAPSD
for short.

Using these two measures, we can analyze how LS performs with different levels
of noise. As can be seen in Fig. 2B, E for both the FitzHugh-Nagumo and
Schnakenberg models and in Fig. S1A and S2A, at relative noise levels below
around 0.5% we obtain very similar patterns, even when the corresponding
relative error in parameters is around 0.1%. For larger levels of noise, both
examples in Fig 2B begin to fail, resulting in parameters that do not produce
patterns. In Fig. 2E we can see that the standard error in parameters (the
shaded region encompassing the solid lines) behaves very differently for the
two models, and its width seems to reduce at different values of relative
error, so it is not caused by the scale of the error. We found that the point
where this error is reduced is when the LS results are not centered around the
original set of parameters anymore. To visualize this, it is helpful to note
that the LS estimator is dependent on noise. Hence, different realizations of
the corrupted pattern at the same level of noise will produce different LS
estimators, which will produce different patterns, as can be seen with $1\%$
noise in Fig. 2B. We can either obtain a pattern or a constant solution,
showing how the obtained parameters are at the boundary of the Turing region.
These different estimators can be thought of as a set of points in parameter
space. Subsequently, the reduction in standard error occurs as the set of
points moves further away from the original set of parameters, as described in
the sketch in Fig. 2E). Note that since we are using LS, instead of computing
the relative error in $\hat{\beta}$, it could be argued that we can use the
statistical properties of the LS estimator and obtain an expression for its
variance. Indeed, the relative error we are computing is proportional to the
standard deviation of $\hat{\beta}$. Hence, if we had an expression for the
variance, then the errors would be easily computed. However, this is not an
easy task since the usual LS assumptions do not hold in our case (noise is not
additive and does not have a constant variance, and the design matrix is
noisy). Applying the statistical properties of the LS estimator results in
errors that do not match our data, so we decided to numerically estimate the
relative error.

In Figs. 2D and S2B we notice discrete jumps when comparing RAPS profiles,
which we can match to behaviors in the patterns. By comparing with different
noise levels in Fig. 2B, we found that the last jump separates inferred
parameters that produce a pattern and inferred parameters that do not. The
other jumps are less clear and have to do with the patterns changing
wavelength and scale. For example, for the Schnakenberg model we can notice a
slow linear trend, which shows that the range of the patterns is slowly
changing. This is followed by a jump around $0.4\%$, which is the point where
the spots have grown too much and the pattern changes to fewer dots (see
difference in spot number between $0.3\%$ and $0.6\%$ noise). For the FitzHugh
Nagumo model, we notice more jumps and it is less clear what these represent.
We saw that the general behavior is conserved for different patterns with the
same parameters, but the position of jumps seem to be dependent on the initial
conditions. Note that the main advantage of this measure is that it allows us
to see when changes in the pattern occur. Nevertheless, this measure is
computationally more expensive than computing the relative error in
parameters, since we need to numerically solve for the pattern. Hence, we
still use the relative error in parameters.

Our results show that the LS estimators are not robust to the levels of noise
we aim for, but in the absence of noise the method works well with little
data. As can be seen in Fig. 3A, even if we crop the original pattern and use
a $4\times 4$ pixel region, we still obtain the same patterns. A natural
question is to find the minimum number of pixels. Given the simplicity of the
method, to answer this question we need only consider the formulation of LS. A
full explanation is provided in the Methods section but we give a brief
explanation here: given the design matrix $X$, which is determined by the
model, we need $X^{T}X$ to be invertible. This condition will be satisfied
when $X$ has enough rows (a pixel corresponds to two rows) so that it is a
full rank matrix, making the minimum amount of pixels be model dependent. For
the Schnakenberg and Brusselator models, 2 pixels are enough, whereas 3 pixels
were required for the FitzHugh-Nagumo model.

Next, we investigated how much the addition of new pixels changes the relative
error in the parameters. To do this, we sampled $N$ pixels from the pattern
using a two dimensional uniform distribution (Fig. 3B), which we used to infer
a new set of parameters. In Fig. 3C we show a scatter plot of the mean
relative errors for each $N.$ We observed that the effect of the number of
pixels on the relative error is different depending on whether the original
pattern was corrupted with noise or not. Without noise (blue) we hardly see
any improvement, but with noise (orange) there is a more drastic improvement
in the parameters, which is expected as the noise has less effect the larger
the sample. As the amount of randomly selected pixels is increased, the
relative error approaches that of no added noise with a power law, as can be
seen in the log-log plot in the upper right corner of Fig. 3C. The exponent of
this power law is close to $-\frac{1}{2}$, which agrees with previous results
on the LS convergence with $N$ [28]. For an intuitive mathematical explanation
of this convergence see Supplementary Information Section 4. We remark that
the fluctuations in relative error without noise (or for large $N$ with noise)
are numerical instabilities, since they are also observed if we use all the
data and merely change the order of the input. Similarly to this, we also
investigated this effect in the RAPSD for the Schnakenberg model in Fig. S3.
We saw that the observed behavior was dependent on the initial level of noise
chosen: for a noise around $0.03\%$, we observe a similar power law as with
the relative error, with a slope again close to $-\frac{1}{2}$. For a noise of
$0.25\%$, we observed the same jumps as can be seen in Fig. 2D, with a
convergence towards the same value shown in the figure for this level of
noise. (For more discussion see Supplementary Information Section 5.)

We also attempted to apply LS to a design matrix defined by a model not used
to generate the pattern. For example, we produce a pattern with the
Brusselator model but we define the LS minimization using the Schnakenberg
model. This is what we called ‘mixing’ the models, and the goal is to check if
the models are flexible enough to produce the same pattern from one another.
We found that this did not work even between models capable of producing the
same type of pattern, and instead produced parameters that gave no Turing
patterns at all (see the Discussion section). We can see that this result
enforces the idea that a model is like a key, and when we have the key, the
pattern is simply a barcode which is easy to read and gives us the model
parameters. In summary, the LS method works well without noise in the the
data, as it returns the exact parameters (with some numerical error), even
when we have very little data (2-3 pixels depending on the model).

Figure 4: Physics-informed neural networks for parameter inference. (A)
Architecture of RBF-PINNs, where the input is space coordinates ($x,y$) and
the output is the pattern. Input is shown in green, variables which are
trained in blue and input to the losses in yellow. Red arrows show denote
usage in the loss and yellow arrows backpropagation. From the network the
partial derivatives can be efficiently computed using automatic
differentiation and used in the PDE loss, where the PDE parameters are also
network parameters. (B) Illustration of the three parameters of the Gaussian
kernel and their interpretation. (C) Results from RBF-PINNs with different
levels of added noise. After adding noise to the Turing pattern, the network
is used to obtain a parameter set, which is subsequently used to predict the
pattern. Patterns to the left are the original ones. (D) Relative error in the
parameters and the RAPS difference for the parameters and patterns used for
(C).

### Physics-informed neural networks

Physics-informed neural networks (PINNs) are a neural network architecture in
which the network serves as a function approximator with two loss functions:
the first compares the output of the network to the numerical simulation of
the PDE (the pattern in our case) and makes sure the network outputs the
correct values; the second tries to make the output of the network a solution
of the PDE (the physical law) by both optimizing the network and optimizing
the parameters of the PDE (the estimation)[25]. Using these two losses, the
network approximates the solution as well as learns the optimal set of
parameters for which the Turing pattern is a steady-state solution. This
approach is computationally very efficient, since it only adds a few
parameters to the network which usually trains thousands. Hence, the extra
computational cost is minimal, leading to excellent scaling with the number of
parameters in the PDE. A peculiarity of this method is that if we have no
noise in our input patterns, we can let the approximation overfit the data.
This is because the approximation will become better and this will improve the
parameter optimization as well.

As function approximation we used a radial basis function neural network
(RBFNN), since it is especially suitable for data with regularities such as
evenly spaced peaks and valleys like our Turing patterns. This network only
has one hidden layer (aside from input and output) defining the kernel. We can
think of each kernel as adding a (unnormalized) distribution (Gaussian in our
case) at a given location with a given variance and weight. These three are
the only types of parameters of this network. Hence, the amount of nodes is
the amount of kernels that we have, and during training these will change
their location and variance, shrinking or growing to approximate the pattern.
A representation of the three parameters and their effects are shown in Fig.
4B.

As with our previous method, without added noise the network can approximate
the pattern up to an arbitrary accuracy [24]. This means that we can reach the
same levels of accuracy as with the other method, but with a caveat: we need
to train the network for a long time to reach this level of accuracy and we
would need a large number of nodes in the hidden layer. This is the reason
behind the large error in Fig. 4D without added noise, since we use the same
amount of nodes for all noise levels. Note the performance can easily be
improved without noise, but as our main focus is robustness to noise, we did
not tune any parameters to any specific level of noise, except for the number
of nodes and the variance of the kernels, and we use exactly the same network.
Once we apply noise, our results with this technique show an improved
robustness when compared to LS, as can be seen in Fig. 4C. For both the
FitzHugh-Nagumo (top) and Schnakenberg (bottom), we obtain very similar
patterns up to a noise level of 12-16%. For the Brusselator model, shown in
Fig. S1B, we can see similar results. Comparison with Fig. 4D and Fig. S2C and
D shows that the relative error in parameters is higher in the Schnakenberg
than in the FitzHugh-Nagumo model, but this error is much lower in the
Brusselator model, indicating that this model is more robust to noise.
Nevertheless, from Fig. 4C we can see that there is not a direct relationship
between relative error and similarity between the patterns (RAPSD), since for
low noise levels (0-4%) the relative error in the Schnakenberg model is higher
and the RAPSD is lower than the FitzHugh-Nagumo model, but for higher levels
the RAPSD of the Schnakenberg model becomes higher. For convergence plots of
some of the model parameters and losses of the network see Fig. S4 in
Supplementary Information Section 6.

By looking at the RAPS profiles we find the same discrete jumps that we saw
with the LS method. In particular, a well pronounced final jump is observed,
which we argued before indicates the transition from pattern to no pattern. By
comparing Figs. 2D and 4D, RBF-PINNs with 0-4% noise results in RAPS values
corresponding to around 0.3-0.5% for LS while the points with 8-12% noise
correspond to the ones around 0.6-0.9%. This means that RBF-PINNs can infer
parameters and predict patterns with around 10-20 times more noise than LS for
the Schnakenberg and FitzHugh-Nagumo models. In the case of the Brusselator
this number is even higher, as can be seen by comparing Fig. S2B and D, with
RBF-PINNs being able to predict patterns with up to 40 times higher noise
levels.

### Application to data from chemical patterns

Having demonstrated how inference by PINNs shows strong robustness to noise,
the next step is to apply this method to experimental data. We chose chemical
Turing patterns as an example because patterns are more stable and robust than
biological patterns, and corresponding models are well established.
Specifically, we consider here the chlorine dioxide–iodine–malonic acid
(CDIMA) system used to study the impact of 2D growth on patterns formation
[8]. This reaction shows a photosensitivity that was utilized to produce a
time series of snapshots of the pattern at different times in a radially
growing domain using a mask. Since our method focuses on the steady state
pattern, we discarded all the time points except for the last one, focusing on
the stable central region away from the boundary.

The experiments were modeled using a Lengyel-Epstein two-variable model [29],
modified to incorporate the effects of illumination [30]. Since we are only
interested in the pattern, which is shown in the dark region, we will use the
original model without illumination. We found that this model, although
already non-dimensional, was problematic for parameter inference. This was
because in one of the equations all terms have a parameter, which can lead to
a trivial solution where all parameters are zero. To solve this, we derived a
new non-dimensionalization:

$u_{t}=d\Delta u+c_{1}-u-c_{2}\frac{4uv}{1+u^{2}}\qquad v_{t}=\Delta
v+u-c_{2}\frac{4uv}{1+u^{2}},\\\ $ (9)

which we used to fit the data with our RBF-PINN. Before stating our results,
there are a few complications worth mentioning. First and foremost, the scale
of the patterns is unknown, since the original data is an image with pixel
values ranging from 23 to 255. Secondly, we have two concentrations on the
model ($u$ and $v$) exhibiting a pattern, but only one experimental pattern.
In order to solve both these problems, we define a free-scale variable $W$,
which is a rescaled version of the original pattern in the range $[0,1]$. In
order to obtain the model concentrations, we assume that there is a linear map
from $W$ to $u$ and $v$, which we can write as:

$u=W\kappa_{u}+\gamma_{u}\qquad v=W\kappa_{v}+\gamma_{v},\\\ $ (10)

where $\kappa_{x}$ and $\gamma_{x}$ for $x=u,v$ are scale and shift parameters
respectively. The assumption of the existence of this linear map can be
justified by the fact that patterns usually are either in phase (positive
$\kappa$) or out of phase (negative $\kappa$).

As before, our network will consist of two independent subnetworks which will
approximate $W$ individually. From these two approximations we will recover
$u$ and $v$ using the scaling in Eq.10, so we will call these approximations
$W_{u}$ and $W_{v}$, respectively. The reason for using two approximations for
the same variable is that usually patterns are not perfectly in phase or out
of phase. Hence, by separating these into two variables we allow for
flexibility in the method to make adjustments to each pattern individually.
Another difference with the network used previously is that the scaling
parameters are added as part of our PDE loss, and hence are optimized with the
rest of parameters, as portrayed in Fig. 5B.

We firstly tried our inference method using numerical ‘data’ from a simulation
of the pattern. This gave us two patterns, but we only used one of them to
make this more similar to the experimental case. As a proof of concept, we
started by fixing the scaling parameters to the best combination possible for
both $u$ and $v$, and obtaining the rest. The predicted pattern is shown in
the first column of Fig. 5C (bottom), together with the original one (top),
and we can see a good correspondence between them. Specifically, both show
labyrinths on a similar scale. We obtained a RAPS of 0.2941 for the two
numerical patterns, which is close to the value obtained for the FitzHugh-
Nagumo model in Fig. 4D below 16% of noise.

Figure 5: Application to chemical patterns. (A) Explanation of scaling
procedure, showing numerical and experimental patterns and how the scaling to
the free-scale variable and the rescaling using the shift and scale parameters
are performed. (B) Architecture of RBF-PINNs for the experimental case, with
the division into the $u$ and $v$ approximation, the rescaling and the
different losses. (C) Results from RBF-PINNs to the numerical pattern (first
column) and the experimental pattern (second column). The top images show the
original patterns, while the bottom images show the predicted patterns using
the inferred parameters and our numerical solver. (D) Time evolution of
simulated experimental pattern showing how labyrinths are present at the
initial time points and plot of convergence of scaling parameters.

Encouraged by the numerical result, we tried our approach on the experimental
pattern. We initialized all parameters close to the values from the numerical
pattern and also used the same spatial dimensions. The results are shown in
the second column of Fig. 5C. We can see that the predicted pattern (bottom)
is not particularly close to the experimentally observed one (top). The
obtained pattern also seems to show spots instead of labyrinths, like the
experimental one. Nevertheless, when inspecting the time evolution of this
pattern in Fig. 5D, we can observe that at the beginning the spots are
connected and later on they move apart. By looking at the experimental videos,
we can actually see this happening too - at the beginning we can only observe
labyrinths but later on spots begin to appear. This seems to point to the fact
that the experimental pattern has not completely converged and that if left to
evolve for longer, it would end up in a steady state with only spots.
Furthermore, in Fig. 5D (bottom) we see the convergence of the scaling
parameters for $u$ (solid lines) along with the values from the original
numerical pattern (dotted lines). We can notice that they match quite well,
with a slight mismatch of the shift parameter which is probably due to
feedback in the training with other parameters. We could also infer that the
chemical pattern that was produced experimentally was most likely an
activator, $I^{-}$, which we later confirmed [31]. This is because depending
on whether the initial concentration belongs to $u$ or $v$ in the model, one
non-dimensionalization will work better than the other; in our case only $u$
worked. The network training is further discussed in the Supplementary
Information section 7.

In summary, our method can be applied to experimental data with only minor
manual intervention. Furthermore, we were able to obtain insight into the
chemical pattern e.g. its closeness to a spotty steady state, due to the
parameters we obtained from the method. Our method predicted the correct
identification of the pattern with the activator chemical species.

### Discussion

We investigated two methods to solve the inverse problem in Turing patterns.
The first was based on the least squares (LS) method, and was found to work
very accurately, with relative errors under $10^{-15}$, when no noise is added
to the patterns. We also proved that it works even with very small quantities
of data: 2-3 pixels from the discrete patterns were enough, and that this
depends on the model equations. We also proposed another method based on
neural networks, that we called RBF-PINNs since it is a mixture of PINNs and
RBF neural networks. Without noise, this method can be compared with LS, since
the universal approximation theorem assures that we are able to approximate
any function up to any given accuracy [24]. The main advantage of this method
comes when we consider a noisy pattern, since in this case our method was
capable of recovering similar parameters to the original ones with up to
12-16% relative noise.This allowed us to use this method for experimental
chemical patterns, and obtain a set of parameters that gave us new insight
about the possible future behavior of the patterns obtained experimentally.

Each of these methods has its strengths and weaknesses. The main disadvantage
of the LS method is its sensitivity to noise, which is usually present in most
biological systems and data. This was remedied by our RBF-PINN method at the
expense of a computational cost several orders of magnitudes larger than LS,
taking nearly an hour to train on a HP Z8 G4 Workstation. Interestingly, the
inferred parameters without noise are not as accurate as LS, which can be
observed by comparing the relative errors in Figs. 2D and 4D. This performance
can be further optimized by using more nodes in the network. Also the
robustness to noise can be improved. Both our approaches used a square grid of
$50\times 50$ pixels, but as we increase the amount of pixels, the performance
of both methods in presence of noise becomes more robust. We can also run the
optimization of the loss functions for longer to obtain a better convergence,
or we can alter the architecture to make it more specific to our patterns. In
this paper, we use exactly the same network architecture for all patterns,
only changing the variance of the starting RBF kernels since the patterns have
a different wavelength and scale, and the number of nodes. However, we could
get more accurate results if we were to increase the number of nodes in the
network, although we would have to be careful about overfitting to noisy data.

When comparing with other attempts at solving the inverse problem for Turing
patterns, we should first make the distinction between approaches that treated
noise and those which did not. For the approaches that did not consider noise
[18], the simple LS performs better than this Bayesian-based approach since it
does not need any training and it is more generalizable, since we only need
knowledge of one pattern and not a library of patterns. Compared to the ones
that did consider noise [21, 19], the main difference is that we do not
require knowledge of the initial conditions or of transient data (dynamics).
These are assumptions that would not hold for biological data, which usually
is scarce and with high levels of noise, so our approach is to the best of our
knowledge the first one that tackles all of these issues (small data and
noise) at the same time. Recently, the inverse problem was addressed with
neural networks but with some limitations: in [32] PINNs were applied to
Turing patterns but without considering noise and in [33] an architecture
based on recurrent neural networks was developed but without testing it on
actual data.

Despite the successes of our methods to deal with small data and noise, there
are some limitations worth mentioning. First, we did not use domain growth for
our work, which would be more biologically relevant. This could be implemented
by considering more traditional PINNs that take time into account. We did not
focus on time series data, although our LS approach easily extends to this
(see Supplementary Information section 9), and the more traditional PINN
structures are very adequate for data changing in time. Our approach to
incorporate the scaling in the chemical patterns, although effective for this
case, was found to be very sensitive, since parameters could interfere with
other parameters during training. A more sophisticated encoding of the scaling
used could potentially provide a better distinction in the parameters and
yield better and more robust results. Looking at possible biological
applications, a problem would be the simplicity of the models, which are just
cartoon representations of the actual biological networks as encountered e.g.
in developmental biology. To inform biological experiments, future models need
to be complex enough to have meaningful connections to experimental
observables, while being simple enough to avoid the ‘curse of dimensionality’.
At the same time, experiments should show a good connection with the
simulations in the models, as otherwise the approaches shown in this work
would not be useful. As a remark, we only focused on finding the parameters,
but it would also be interesting to extend this to model selection. Existing
methods include graph networks [34], embedding theory and diffusion maps [35,
36].

Our approaches open up new ways of connecting mathematical models to
experimental patterns. In particular, RBF-PINNs worked for noise levels
comparable to biological data, and proved to be useful at elucidating
properties of chemical patterns. Hence, this method could potentially be
applied to infer parameters for biological candidate models, “proving” that a
model is capable of reproducing observed patterns. This type of model
selection could ultimately aid the rational design of synthetic tissues with
patterns for downstream templating and added functionality [37, 38]. In
conclusion, we hope our machine learning approaches to solving inverse
problems stimulate new research into unraveling pattern formation in biology,
chemistry, and bioengineering.

## Methods

All algorithms and methods were run in Python on a HP Z8 G4 workstation with a
Intel® Xeon(R) Gold 6128 CPU @ 3.40GHz × 24 processor and a Quadro RTX
6000/PCIe/SSE2 GPU. Packages and requirements of the environments are
explained in the Supplementary Information section 8.

### Numerical simulations

All numerical simulations of the models were performed on a $50\times 50$
discrete grid by applying a center difference approximation to the second-
order derivatives in the Laplacian, hence transforming the PDE Eq. 1 into an
ODE model:

$\displaystyle\begin{split}\frac{du}{dt}=&D_{u}\frac{u_{i,j-1}+u_{i,j+1}+u_{i-1,j}+u_{i+1,j}-4u_{i,j}}{(\Delta
x)^{2}}+f(u,v)\\\
\frac{dv}{dt}=&D_{v}\frac{v_{i,j-1}+v_{i,j+1}+v_{i-1,j}+v_{i+1,j}-4v_{i,j}}{(\Delta
x)^{2}}+g(u,v)\\\ \end{split}$ (11)

Here, we assumed that the step in the $x$ direction is the same as the step in
the $y$ direction, i.e $\Delta x=\Delta y$. With this approximation, we use a
numerical solver to forward integrate these equations for a sufficiently large
amount of time that ensures that the system has converged to a steady state,
which is the Turing pattern.

### Least squares

As described in the section on LS, in order to write down the solution we
first need to write our problem in matrix form. We will do this here with the
Schnakenberg model as an example. Using the same notation as before, we let
$u_{ij}$ and $v_{ij}$ be our discretized Turing patterns in matrix form, for
$i,j=1,2...,N$, with $N$ being the number of columns and rows (as we have a
square grid). Writing Eq. 5 using this notation and at steady state we obtain:

$u_{ij}-\Delta u_{ij}=c_{1}+c_{2}u_{ij}^{2}v_{ij}\qquad 0=d\Delta
v_{ij}+c_{3}-c_{2}u_{ij}^{2}v,\\\ $ (12)

where we collected the terms not multiplied by a parameter to the left. This
is a total of $2\times N^{2}$ equations, which we can write in matrix form as:

$\begin{pmatrix}u_{11}-\Delta u_{11}\\\ u_{12}-\Delta u_{12}\\\ \vdots\\\
u_{NN}-\Delta u_{NN}\\\ 0\\\ 0\\\ \vdots\\\
0\par\end{pmatrix}=\begin{bmatrix}1&u_{11}^{2}v_{11}&0&0\\\
1&u_{12}^{2}v_{12}&0&0\\\ \vdots&\vdots&\vdots&\vdots\\\
1&u_{NN}^{2}v_{NN}&0&0\\\ 0&-u_{11}^{2}v_{11}&1&\Delta v_{11}\\\
0&-u_{12}^{2}v_{12}&1&\Delta v_{12}\\\ \vdots&\vdots&\vdots&\vdots\\\
0&-u_{NN}^{2}v_{NN}&1&\Delta v_{NN}\par\end{bmatrix}\begin{pmatrix}c_{1}\\\
c_{2}\\\ c_{3}\\\ d\end{pmatrix}$

We will write this equation as $y=X\beta$, and will call $y$ the vector of
independent variable and $X$ the matrix of dependent variables. It is also
customary when performing regression to call $y$ the vector of outputs and $X$
the design matrix. This is the starting form of the LS formulation [22]. By
defining the error vector and finding its minimum, the solution for the
optimal parameters can be written as:

$\beta=(X^{T}X)^{-1}X^{T}y$

Note that for this approach to work we need $X^{T}X$ to be invertible, which
is equivalent to $X$ having full rank. Based on this necessary condition, we
can obtain the minimum amount of pixels needed for the solution to be well
defined. In the case of the Schnakenberg model, we can see that two pixels are
enough, since it gives us a $4\times 4$ matrix in the form:

$M=\begin{bmatrix}1&a&0&0\\\ 1&b&0&0\\\ 0&-a&1&c\\\ 0&-b&1&d\end{bmatrix}$

which for most $a,b,c,d$ will be full rank. There are cases (e.g. $c=d$ or
$a=b$) where $M$ will not be full rank, but for randomly selected pixels in
the pattern these cases are unlikely. Likewise, we can easily see that this
condition gives two pixels for the Brusselator model (since we have two
independent $2\times 2$ submatrices) and three for the FitzHugh-Nagumo model
(since we have a submatrix of $3\times 3$).

### Neural networks

For the neural network, as described in the Results section, we combined PINNs
(i.e. using the PDE equations in order to improve the approximation and having
the parameters as trainable variables in the network) with RBF neural networks
to obtain what we called RBF-PINNs. The architecture is very simple since the
network is made of only three layers. The input of the network is a vector
representing a spatial location, in our case $(x,y)$, since we only have two
dimensions, and the output is $u(x,y)$. We have one network for each $u$ and
$v$, to allow for more flexibility in the model, but we train them together.
Between the input and the output layer we have a single hidden layer, with an
activation layer representing the RBF kernels, and the input being fed
directly to this activation layer by having all weights set to 1, which means
that the amount of trainable parameters is in the order of $3N$ if we have $N$
nodes. These RBF kernels can be written as:

$\displaystyle\begin{split}\phi_{i}(\mathbf{x})=w_{i}e^{-\beta_{i}||\mathbf{x}-\mathbf{c_{i}}||^{2}}\quad\text{for}\quad
i=1,\>2,...\>M\end{split}$ (13)

where $M$ is the total number of nodes in the network. As explained previously
and portrayed in Fig. 4, $\beta_{i}$ works as the variance, $c_{i}$ is the
location of the $i^{th}$ kernel and $w_{i}$ its weight or importance. Each
node in the network represents a different kernel, so the only way to change
the architecture of the network to improve the accuracy of the approximation
is to increase the amount of kernels by increasing the amount of nodes. The
kernel parameters were set randomly, $\mathbf{c}$ was set to a 2D uniform
random variable in the spatial range of the pattern, and the network weights
were initialized using Glorot initialization. The variance was initialized
depending on the pattern, since the variance of the Gaussian kernels has to be
on a similar scale as the pattern wavelength.

Our network has two main losses. One is an approximation loss, which compares
the output of the network to the pattern. If we consider our network as a
function, we can write it as $\Phi^{u}(\textbf{X})$ where $\mathbf{X}$ is the
space dimensions where we have the pattern values $\mathbf{u}$. We will from
now on just write this as $\Phi^{u}$ and similarly for $\mathbf{v}$.
Subsequently, we can write the approximation loss for $\mathbf{u}$ as:

$\displaystyle\begin{split}L_{app,u}=\frac{1}{N^{2}}\sum_{i,j}||\Phi^{u}_{i,j}-\mathbf{u}_{i,j}||^{2}\end{split}$
(14)

Note that this is equivalent to computing the mean squared error (MSE) between
the data and the approximation. We mentioned the RAPS loss before, and it
could be used here for the network training as well, but since we care about
the pattern being accurate, the MSE works fine. However this could be a
possible improvement since it should cancel the noise to some extent.

The second loss is the PDE loss, which enforces that the output has a small
PDE residual and at the same time improves the parameters so that the data
gives a better fit to the equation:

$\displaystyle\begin{split}L_{PDE}=\frac{1}{N^{2}}\sum_{i,j}||\beta_{1}\Delta\Phi^{u}_{i,j}+f(\Phi^{u}_{i,j},\Phi^{v}_{i,j},\boldsymbol{\beta})||^{2}+||\Delta\Phi^{v}_{i,j}+g(\Phi^{u}_{i,j},\Phi^{v}_{i,j},\boldsymbol{\beta})||^{2}\end{split}$
(15)

We also have a loss for the diffusion term, to make sure that the second
derivatives of the approximation are accurate. The final structure used for
the results in the paper has 60-120 nodes depending on the model. For
training, we use batch training with batches of 128 elements for 200,000
iterations and Adam optimizer, taking a total of less than an hour to run. The
first 10,000-20,000 iterations are used to approximate the function, since we
only update the weights using $L_{app}$, and after this we change to
minimizing both losses, taking care of using weights before each loss so that
the approximation does not get completely changed by the addition of the new
losses. When training the approximation, we use information on the whole
pattern, but when we switch to the $L_{PDE}$, we only use interior points so
that inaccuracies in the diffusion do not present a problem.

The structure of our network is relatively simple and the RBF network
guarantees a smooth output. However, to avoid overfitting we can use a small
number of nodes or kernels. This is because overfitting is only a problem for
the approximation part of our network (that is where our data is used).
Furthermore, our patterns have a larger wavelength than the noise. By
selecting a variance or amplitude parameter to be of the same order as our
pattern, we keep the network from incorporating most of the noise. On top of
that, we redefine the weights of the network every 2,000 iterations so that
all of them are on a similar scale. This is useful since it is possible that
some weights converge faster than other, and so would eventually have a very
small loss if the original weights were used.

## Acknowledgments

We thank Roozbeh H. Pazuki for invaluable support and technical comments,
Martina Oliver Huidobro for stimulating discussions and Milos Dolnik for
explanations of his experiments. This research was funded through a
studentship from the Department of Life Sciences at Imperial College London.

## Declaration of interests

The authors declare no competing interests.

## Author contributions

A.M.G. and R.G.E. designed, and A.M.G. implemented and performed the
theoretical computational approach and data analysis. Both authors wrote the
paper.

## References

  * [1] Turing AM.  The Chemical Basis of Morphogenesis.  Philosophical Transactions of the Royal Society of London Series B, Biological Sciences. 1952;237(641):37-72.  Available from: http://www.jstor.org/stable/92463.
  * [2] Schnakenberg J.  Simple chemical reaction systems with limit cycle behaviour.  Journal of Theoretical Biology. 1979 12;81:389-400.
  * [3] Gierer A, Meinhardt H.  Gierer, A. and Meinhardt, H. A theory of biological pattern formation. Kybernetik 12, 30- 39.  Biological Cybernetics. 1972 01;12:30-9.
  * [4] Prigogine I, Lefever R.  Symmetry Breaking Instabilities in Dissipative Systems. II.  The Journal of Chemical Physics. 1968 2;48:1695-700.  Doi: 10.1063/1.1668896.  Available from: https://doi.org/10.1063/1.1668896.
  * [5] FitzHugh R.  Impulses and Physiological States in Theoretical Models of Nerve Membrane.  Biophysical Journal. 1961;1:445.  Available from: /pmc/articles/PMC1366333/?report=abstracthttps://www.ncbi.nlm.nih.gov/pmc/articles/PMC1366333/.
  * [6] Nagumo J, Arimoto S, Yoshizawa S.  An Active Pulse Transmission Line Simulating Nerve Axon.  Proceedings of the IRE. 1962;50(10):2061-70.
  * [7] Castets V, Dulos E, Boissonade J, De Kepper P.  Experimental evidence of a sustained standing Turing-type nonequilibrium chemical pattern.  Phys Rev Lett. 1990 Jun;64:2953-6.  Available from: https://link.aps.org/doi/10.1103/PhysRevLett.64.2953.
  * [8] Konow C, Somberg NH, Chavez J, Epstein IR, Dolnik M.  Turing patterns on radially growing domains: experiments and simulations.  Phys Chem Chem Phys. 2019;21:6718-24.  Available from: http://dx.doi.org/10.1039/C8CP07797E.
  * [9] Sheth R, Marcon L, Bastida MF, Junco M, Quintana L, Dahn R, et al.  ¡i¿Hox¡/i¿ Genes Regulate Digit Patterning by Controlling the Wavelength of a Turing-Type Mechanism.  Science. 2012;338(6113):1476-80.  Available from: https://www.science.org/doi/abs/10.1126/science.1226804.
  * [10] Watanabe M, Kondo S.  Is pigment patterning in fish skin determined by the Turing mechanism?  Trends in Genetics. 2015;31(2):88-96.  Available from: https://www.sciencedirect.com/science/article/pii/S0168952514001978.
  * [11] Glover JD, Wells KL, Matthäus F, Painter KJ, Ho W, Riddell J, et al.  Hierarchical patterning modes orchestrate hair follicle morphogenesis.  PLoS Biol. 2017 Jul;15(7):e2002117.
  * [12] Marcon L, Sharpe J.  Turing patterns in development: what about the horse part?  Curr Opin Genet Dev. 2012 Dec;22(6):578-84.
  * [13] Painter KJ, Ptashnyk M, Headon DJ.  Systems for intricate patterning of the vertebrate anatomy.  Philos Trans A Math Phys Eng Sci. 2021 Nov;379(2213):20200270.
  * [14] Karig D, Martini KM, Lu T, DeLateur NA, Goldenfeld N, Weiss R.  Stochastic Turing patterns in a synthetic bacterial population.  Proceedings of the National Academy of Sciences. 2018;115(26):6572-7.  Available from: https://www.pnas.org/doi/abs/10.1073/pnas.1720770115.
  * [15] Sekine R, Shibata T, Ebisuya M.  Synthetic mammalian pattern formation driven by differential diffusivity of Nodal and Lefty.  Nature Communications. 2018 Dec;9(1):5456.  Available from: https://doi.org/10.1038/s41467-018-07847-x.
  * [16] Scholes NS, Schnoerr D, Isalan M, Stumpf MPH.  A Comprehensive Network Atlas Reveals That Turing Patterns Are Common but Not Robust.  Cell Systems. 2019;9(3):243-57.e4.  Available from: https://www.sciencedirect.com/science/article/pii/S2405471219302649.
  * [17] Murray JD.  Chapter 1.  In: Spatial Pattern Formation with Reaction Diffusion Systems. New York, NY: Springer New York; 2003. p. 71-140.  Available from: https://doi.org/10.1007/0-387-22438-6_2.
  * [18] Kazarnikov A, Haario H.  Statistical approach for parameter identification by Turing patterns.  Journal of Theoretical Biology. 2020 9;501:110319.
  * [19] Campillo-Funollet E, Venkataraman C, Madzvamuse A.  Bayesian Parameter Identification for Turing Systems on Stationary and Evolving Domains.  Bulletin of Mathematical Biology. 2019 1;81:81-104.  Available from: https://link.springer.com/article/10.1007/s11538-018-0518-z.
  * [20] Schnörr D, Schnörr C.  Learning System Parameters from Turing Patterns. 2021 8.  Available from: https://arxiv.org/abs/2108.08542v1.
  * [21] Garvie MR, Maini PK, Trenchea C.  An efficient and robust numerical algorithm for estimating parameters in Turing systems.  Journal of Computational Physics. 2010;229(19):7058-71.  Available from: https://www.sciencedirect.com/science/article/pii/S0021999110003128.
  * [22] Hastie T, Tibshirani R, Friedman J.  The Elements of Statistical Learning.  Springer Series in Statistics. New York, NY, USA: Springer New York Inc.; 2001.
  * [23] Hornik K, Stinchcombe M, White H.  Multilayer feedforward networks are universal approximators.  Neural Networks. 1989;2(5):359-66.  Available from: https://www.sciencedirect.com/science/article/pii/0893608089900208.
  * [24] Park J, Sandberg IW.  Universal Approximation Using Radial-Basis-Function Networks.  Neural Computation. 1991;3(2):246-57.
  * [25] Raissi M, Perdikaris P, Karniadakis GE.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics. 2019;378:686-707.  Available from: https://www.sciencedirect.com/science/article/pii/S0021999118307125.
  * [26] Cavanagh H, Mosbach A, Scalliet G, Lind R, Endres RG.  Physics-informed deep learning characterizes morphodynamics of Asian soybean rust disease.  Nature Communications. 2021 Nov;12(1):6424.  Available from: https://doi.org/10.1038/s41467-021-26577-1.
  * [27] Maini PK, Woolley TE, Baker RE, Gaffney EA, Lee SS.  Turing’s model for biological pattern formation and the robustness problem.  Interface Focus. 2012 8;2:487-96.  Available from: https://royalsocietypublishing.org/doi/full/10.1098/rsfs.2011.0113.
  * [28] van de Geer S.  Estimating a Regression Function.  The Annals of Statistics. 1990;18(2):907-24.  Available from: http://www.jstor.org/stable/2242140.
  * [29] Lengyel I, Epstein IR.  Modeling of Turing Structures in the Chlorite-Iodide-Malonic Acid-Starch Reaction System.  Science. 1991;251(4994):650-2.  Available from: https://www.science.org/doi/abs/10.1126/science.251.4994.650.
  * [30] Muñuzuri AP, Dolnik M, Zhabotinsky AM, Epstein IR.  Control of the Chlorine Dioxid-Iodine-Malonic Acid Oscillating Reaction by Illumination.  Journal of the American Chemical Society. 1999;121(35):8065-9.  Available from: https://doi.org/10.1021/ja9910457.
  * [31] Preska Steinberg A, Epstein IR, Dolnik M.  Target Turing patterns and growth dynamics in the chlorine dioxide-iodine-malonic acid reaction.  J Phys Chem A. 2014 Apr;118(13):2393-400.
  * [32] Kho J, Koh W, Wong JC, Chiu PH, Ooi CC. Design of Turing Systems with Physics-Informed Neural Networks; 2022.
  * [33] Rao C, Ren P, Wang Q, Buyukozturk O, Sun H, Liu Y.  Encoding physics to learn reaction–diffusion processes.  Nature Machine Intelligence. 2023 Jul;5(7):765-79.  Available from: https://doi.org/10.1038/s42256-023-00685-7.
  * [34] Cranmer M, Sanchez-Gonzalez A, Battaglia P, Xu R, Cranmer K, Spergel D, et al.. Discovering Symbolic Models from Deep Learning with Inductive Biases. arXiv; 2020.  Available from: https://arxiv.org/abs/2006.11287.
  * [35] Lee S, Kooshkbaghi M, Spiliotis K, Siettos CI, Kevrekidis IG.  Coarse-scale PDEs from fine-scale observations via machine learning.  Chaos: An Interdisciplinary Journal of Nonlinear Science. 2020;30(1):013141.  Available from: https://doi.org/10.1063/1.5126869.
  * [36] Lee S, Psarellis YM, Siettos CI, Kevrekidis IG. Learning black- and gray-box chemotactic PDEs/closures from agent based Monte Carlo simulation data. arXiv; 2022.  Available from: https://arxiv.org/abs/2205.13545.
  * [37] Davies JA, Cachat E.  Synthetic biology meets tissue engineering.  Biochemical Society Transactions. 2016 Jun;44(3):696-701.  Available from: https://doi.org/10.1042/bst20150289.
  * [38] Toda S.  Synthetic tissue engineering: Programming multicellular self-organization by designing customized cell-cell communication.  Biophysics and Physicobiology. 2020;17(0):42-50.  Available from: https://doi.org/10.2142/biophysico.bsj-2020002.


# Forecasting Post-Wildfire Vegetation Recovery in California using a
Convolutional Long Short-Term Memory Tensor Regression Network

Jiahe Liu
Western Connecticut State University
Danbury, CT 06810
<EMAIL_ADDRESS>
Xiaodi Wang
Western Connecticut State University
Danbury, CT 06810
<EMAIL_ADDRESS>

###### Abstract

The study of post-wildfire plant regrowth is essential for developing
successful ecosystem recovery strategies. Prior research mainly examines key
ecological and biogeographical factors influencing post-fire succession. This
research proposes a novel approach for predicting and analyzing post-fire
plant recovery. We develop a Convolutional Long Short-Term Memory Tensor
Regression (ConvLSTMTR) network that predicts future Normalized Difference
Vegetation Index (NDVI) based on short-term plant growth data after fire
containment. The model is trained and tested on 104 major California wildfires
occurring between 2013 and 2020, each with burn areas exceeding 3000 acres.
The integration of ConvLSTM with tensor regression enables the calculation of
an overall logistic growth rate $k$ using predicted NDVI. Overall, our
$k$-value predictions demonstrate impressive performance, with 50% of
predictions exhibiting an absolute error of 0.12 or less, and 75% having an
error of 0.24 or less. Finally, we employ Uniform Manifold Approximation and
Projection (UMAP) and KNN clustering to identify recovery trends, offering
insights into regions with varying rates of recovery. This study pioneers the
combined use of tensor regression and ConvLSTM, and introduces the application
of UMAP for clustering similar wildfires. This advances predictive ecological
modeling and could inform future post-fire vegetation management strategies.

##  1 Introduction

Wildfires are one of the most damaging and costly natural disasters in the
United States. Data released by the National Interagency Fire Center (NIFC)
indicates that there has been a gradual rise in the moving averages of both
total burned acreage and burned acreage per fire over the past few decades
[1]. From 2012 to 2022, an average of 7.3 million acres of vegetation were
burned by wildfires annually in the U.S., with 3 extreme years from 2015 to
2021 where the yearly burned acreage exceeded 10 million acres, surpassing all
pre-2015 records in modern data collection [1]. This results in heightened
costs, ranging from direct costs such as firefighting efforts and property
damage, to indirect costs such as the resulting implications for water
supplies and public health. In the U.S., dedicated federal funding for fire
suppression alone reached $3.65 billion as of 2020 [2].

Across the US, California is by far impacted the most by wildfires, which has
had more than 59,000 structures destroyed by wildfires between 2005 and 2020
[2]. Out of the 20 costliest and most destructive fires in California’s
history, 10 of them have taken place since 2015 [3]. Because California is a
hotspot for major wildfires, it serves as an area of interest for our study.

Post-fire succession refers to the natural sequence of plant and ecosystem
changes that occur after a wildfire, including the colonization of plant
species, regrowth of vegetation, and reestablishment of a

Figure 1: A scatter plot of the coordinates of every wildfire collected in
this study. The large wildfires represented by red dots are included in our
dataset, while the small wildfires represented by blue dots are not.

stable ecosystem. Human intervention may be required to hasten the post-fire
succession of severely damaged forests, as it can be difficult for these areas
to naturally recover [3]. Therefore, the study of post-fire plant regeneration
plays a crucial role in guiding decision-making in conservation efforts.

In this study, a Convolutional Long Short-Term Memory (ConvLSTM) network was
utilized to forecast future vegetation levels up to two years after the
containment of large wildfires in California. Then, we utilized logistic curve
fitting and tensor regression to derive logistic curve parameters, namely a
growth rate $k$ for each fire, creating a combined ConvLSTM Tensor Regression
(ConvLSTMTR) network. Comparing the fitted $k$-values on actual post-wildfire
vegetation levels with the fitted $k$-values on our ConvLSTMTR results, we
found that the framework performed consistently well on out-of-sample
wildfires, with 50% of predictions having an absolute error of 0.12 or less.

###  1.1 Application Context

Many studies around post-wildfire forest regeneration have been conducted by
analyzing specific species of trees or vegetation, allowing us to understand
it from an ecological and biogeographical perspective. The severity of a
wildfire can significantly damage seed availability [4, 5] and forest soils,
slowing the rate of recovery [6]. These factors all impact the duration of
post-wildfire regeneration, which can take years to fully occur [8, 9]. Deep
learning has been explored in a study that forecasts fire danger maps [10].
However, knowing how to respond to the aftermath of a wildfire is just as
important, especially as they have been increasing in severity and burn area.
We propose a framework for long-term post-wildfire recovery analysis, which
will provide the projections needed for disaster response communities to know
which regions of vegetation need the most resources for recovery.

##  2 Data and Methods

###  2.1 Study Area

We selected the state of California as our study area because of its frequent
and major wildfires, generating a comprehensive dataset that enhances the
depth of our analysis. Data points including coordinates, burn areas, and fire
start and end dates for wildfires from 2013 were publicly available and
retrieved from the California Department of Forestry and Fire Protection (CAL
FIRE). Our goal is to analyze large wildfires whose burn areas were
significant enough so that clear vegetation growth patterns could be analyzed,
so we excluded wildfires with a total burn area of less than $3,000$ acres
(12.14 km2), leaving our dataset with 104 wildfires. Figure 1 shows the
distribution of all wildfires collected and highlights the large wildfires
selected in this study.

###  2.2 Data Collection

In this study, global spatial indices were collected from NASA’s Moderate
Resolution Imaging Spectroradiometer (MODIS) to represent key features of
post-wildfire recovery: vegetation levels, surface temperature, precipitation,
and a mask of vegetation area affected by each wildfire. The target variable
to be predicted was the Normalized Difference Vegetation Index (NDVI)
collected from MOD13Q1 [11]. NDVI is a commonly used index to quantify
vegetation greenness and biomass density that is derived from Landsat surface
reflectance indices. As a supplement to NDVI, the Enhanced Vegetation Index
(EVI) was also collected from MOD13Q1, as it exhibits greater sensitivity to
areas with high biomass than NDVI. This gave our model a wider and more
accurate range of vegetation values. Both NDVI and EVI have a range of $-0.2$
to $1$.

Since temperature plays a major role in plant growth [12, 7], the Land Surface
Temperature (LST) index was collected from MOD11A2 [13]. A fire mask was also
collected from MOD14A2 to indicate the severity of wildfires in affected
areas. Lastly, daily precipitation data was acquired from NASA’s Integrated
Multi-satellite Retrievals for GPM (IMERG) [14]. The length of each time
series was set to 25 time steps with an interval of 32-day interval, starting
from the month when each wildfire was marked as contained. For each of these
indices, a 50 by 50-pixel bounding box was extracted with the coordinates of
each wildfire’s origin as the center.

###  2.3 Data Preprocessing

##### NDVI Value Imputation

As MOD13Q1 derives its NDVI calculations from satellite-detected reflectance
signals, the quality of its data is significantly impacted by fluctuations in
snow and ice cover, cloud shadows, and different viewing angles. This would
introduce noise in our dataset, disrupting our model’s ability to accurately
capture temporal trends and spatial patterns in post-fire vegetation recovery.
To address this issue, we extracted a pixel reliability channel from MOD13Q1
to assign a quality rank to each pixel. Pixel reliability values of 2 or 3
indicate that these regions were affected by snow, ice, or cloud cover, so we
assigned NaN values to replace unreliable pixels. We then utilized KNN
imputation to fill in missing values, allowing us to preserve variability by
taking neighbors from both spatial and temporal dimensions.

##### NDVI De-Seasonalization and Normalization

NDVI values exhibit noticeable seasonality trends due to the natural cycles
and processes that affect plant growth. These trends can differ in both
magnitude and time frame, resulting in the need for normalization across all
NDVI frames collected. We collected a series of monthly reference NDVI frames
up to one year before the start of each wildfire. Then, for every post-
wildfire NDVI frame, its pixel values would be divided by the corresponding
values from the reference NDVI frame of the same month. This approach
transformed our model’s task from forecasting absolute and fluctuating NDVI
values to an NDVI ratio, where a value of 1 signifies pre-wildfire vegetation
levels.

###  2.4 Convolutional LSTM Model Architecture

Figure 2: A comparison between the actual post-fire NDVI frames (top) and our
model’s forecasted NDVI frames (bottom) for the Saddle Ridge Fire of 2019.

The model of choice for this study was a Convolutional Long Short-Term Memory
Network (ConvLSTM). The ConvLSTM model [15] is a neural network architecture
that combines Convolutional Neural Networks [16] with Long Short-Term Memory
[17] networks. It is designed to handle sequences of data that have both
spatial and temporal dependencies, making it particularly useful for
spatiotemporal data processing and time-series forecasting.

The convolutional windows of the model are used to recognize shapes and
patterns in areas of low vegetation levels post-wildfire. The hidden states
and memory of the model are used to learn temporal patterns in vegetation
recovery. The ConvLSTM is uniquely designed to handle spatiotemporal
information because all inputs $x_{1},...,x_{t}$, cell outputs
$c_{1},...,c_{t}$, hidden states $h_{1},...,h_{t}$, and gates
$i_{t},f_{t},o_{t}$ are 3D tensors where the last two dimensions of each
tensor represent its spatial dimensions. Model inputs represent series of
images that contain the variables used in model training, while model outputs
represent series of images that contain the model’s forecasted NDVI values.

In our model architecture, the layer responsible for spatiotemporal learning
was the ConvLSTM2D layer. Three such layers were used with filter counts of
$32,128,$ and $64$, and a 3D convolutional layer with a single filter was
utilized to produce the final output, a time series of forecasted NDVI ratios.
A kernel size of ($3,3$) was selected for each ConvLSTM2D layer, and a kernel
size of ($3,3,3$) was selected for the Conv3D layer. The model utilized a
zero-padding technique and a convolutional stride size of $(1,1)$ so that the
size of the filtered spatial data was identical to that of the input data. We
chose the Rectified Linear Unit (ReLU) function as the activation function for
every convolutional layer. Following the first two ConvLSTM2D layers, each
tensor was normalized using batch normalization layers, which transform the
data along each channel to have a mean output value close to 0 and an output
standard deviation close to 1.

##  3 Model Training and Results

###  3.1 Training Data

The input to our ConvLSTM consisted of 5 dimensions, namely samples,
timesteps, rows, columns, and channels. Many of the wildfires included in our
dataset had short-term re-burns either before or after the post-wildfire
analysis period, as well as extreme cases of snow, ice, or cloud cover that
disrupted all or most of some timestamps. These wildfires had NDVI ratio
values that were extremely erratic and unrepresentative of typical a logistic
growth rate, sometimes jumping from average ratio values of 0.5 to 3 or higher
in a single timestamp. Such wildfires would be detrimental to model learning,
and so we selected $31$ wildfires whose average NDVI ratios showed no erratic
increases or decreases to be included in our model’s training and validation.

Given the limited number of wildfires, we then partitioned each original 50 x
50 image into 25 sub-images each spanning 10 x 10 pixels. This increased our
sample size and reduced the complexity of our learning objective by
significantly reducing the size of each sample. Our model was subsequently
able to prioritize the accurate prediction of the NDVI ratios for each pixel,
rather than identifying the larger shapes of the burned areas.

With our dataset, we trained our model for 100 epochs, employing 80% of our
filtered data for training and the remaining for validation. We use a learning
rate scheduler to change the model’s learning rate during training, and
selected mean absolute error as our training and validation loss function.

###  3.2 Results

Once the ConvLSTM model was trained, it was tested by running it on randomly
selected sample wildfires. The model was provided with 5 consecutive time
steps of post-fire observations as input and generated forecasts for the
subsequent 20 frames of NDVI data, covering a recovery period of up to 2 years
after the fire cessation. Figure 2 shows a comparison between the actual and
forecasted NDVI values for the Saddle Ridge Fire of 2019. Our ConvLSTM model
adeptly captures the overall plant vegetation levels, forecasting the temporal
and spatial patterns associated with the data.

##  4 Analysis and Discussion

###  4.1 Plant Growth Curves

While forecasted NDVI time series provide future vegetation index recovery
images, individual NDVI pixels are small and sensitive to random changes and
errors in measuring tools. In order to create a general estimator for
recovery, tensor regression was chosen to generalize the framework, providing
quantified metrics for subsequent analysis and modeling purposes. We chose to
use tensor regression to fit the patterns of our post-wildfire recovery to a
logistic curve, which is commonly used to describe the growth of plants [18].
Similar to a logistic curve, the rate of vegetation recovery is the fastest
immediately after wildfire containment [9]. We first fitted actual post-
wildfire NDVI ratios to gain a baseline set of logistic curve parameters for
each pixel of each wildfire, namely a growth capacity $L$ and a growth rate
$k$. We could then average these values over each wildfire study area and
compare those averages with the logistic growth parameters obtained from our
fitted model predictions.

###  4.2 Convolutional LSTM Tensor Regression Network

Upon obtaining these logistic curve parameters, we enhanced our ConvLSTM model
by incorporating Tucker tensor regression [19] to efficiently translate NDVI
forecasts into measurable logistic growth recovery measurements. This would
extend our pixel-wise logistic curve analysis to a higher-dimensional data
structure, and our model framework was combined to create a Convolutional LSTM
Tensor Regression (ConvLSTMTR) network. Since tensor regression requires
fitting against a collection of tensors along with their target variables’
ground truth values, we used the actual $25$

Figure 3: A histogram of out-of-sample absolute differences between k-values
of logistic fits on actual NDVI data vs. predicted data.

frames of NDVI for each subgrid that had a burn area of at least 50%. We
fitted a logistic curve to every pixel in these grids and averaged the fitted
values for each grid. With this input, the Tucker tensor regression model
could generate a series of estimated logistic curve parameters, providing the
final predictions of our ConvLSTMTR.

We analyzed the 73 out-of-sample wildfires and measured the absolute
difference between the estimated logistic parameters from tensor regression
and those obtained through direct logistic curve fitting on actual NDVI time
series data, as shown in Figure 3. Over 50% of predictions had an absolute
error of less than $0.12$, 75% had an error of less than $0.24$, and 90% of
predictions had an error of less than $0.48$. Outliers occurred in less than
10% of scenarios. These results showed great promise and demonstrated that our
ConvLSTMTR network was not only able to correctly recognize the growth
patterns but was also able to accurately adapt to various rates of vegetation
recovery.

###  4.3 Significance

Previous research on post-fire plant regeneration has primarily focused on
studying plant growth patterns from an ecological and biogeographical
perspective, investigating the mechanisms of post-fire plant recovery, and
analyzing the impact of various factors on this process. However, there has
been a noticeable gap in the literature when it comes to utilizing machine
learning models for spatial-temporal forecasting of plant regrowth
trajectories. Our research addresses this gap by demonstrating the
effectiveness of the ConvLSTM model in accurately predicting future spatial-
temporal plant growth patterns, leveraging various global satellite data.

Considering the vulnerability of machine learning models to data quality
issues, we integrated KNN imputation to replace bad data. Moreover, we
pioneered the application of tensor regression in conjunction with the
ConvLSTM, known as the ConvLSTMTR network, to estimate post-wildfire
vegetation recovery rates. Our approach combines the strengths of the ConvLSTM
in handling high-dimensional spatial-temporal datasets and generating accurate
forecasts with the capabilities of tensor regression to extract growth
features for interpretation and analysis. As a result, our ConvLSTMTR network
exhibits high accuracy in forecasting post-fire recovery rates.

##  5 Conclusion

In this study, we developed a Convolutional Long Short-Term Memory Tensor
Regression network (ConvLSTMTR) to predict plant regeneration in areas
impacted by wildfires. Approximately 100 wildfires in California from 2013
through 2020 with a burn area of at least 3,000 acres were analyzed, and 26
months of post-fire plant recovery data were used for each wildfire. This data
was represented as multi-dimensional spatial-temporal inputs encompassing
variables such as NDVI, EVI, a burned area fire mask, precipitation, and
surface temperature. To effectively address challenges posed by data quality,
limited sample size, and seasonal patterns, we implemented preprocessing
techniques including KNN imputation, stratified sampling, and de-
seasonalization. Upon completing the training stage, the model was able to
forecast NDVI time series 18 months into the future given 5 to 6 months of
immediate post-wildfire observations. Following our ConvLSTM forecasting, we
used Tucker tensor regression to fit logistic growth and capacity parameters
for each predicted NDVI time series. The results from out-of-sample testing
indicated that our ConvLSTMTR was able to forecast a logistic growth rate $k$
for each wildfire with impressive performance. Over 50% of predictions
exhibited an absolute error of less than 0.12, and over 75% had an error of
less than 0.24.

## References

[1] National Interagency Fire Center (2022) Total wildfires and acres.
https://www.nifc.gov/sites/default/files/document-media/TotalFires.pdf.

[2] Troy, A. & Pusina, T. & Romsos, S. & Modhaddas, J. & Buchholz, T. (2022)
The true cost of wildfire in the western U.S. Western Forestry Leadership
Coalition, https://www.thewflc.org/sites/default/files/TrueCostofWildfire.pdf.

[3] California Department of Fish and Wildlife (2020) Science: wildfire
impacts. https://wildlife.ca.gov/Science-Institute/Wildfire-Impacts.

[4] Shive, K. L. & Preisler, H.K. & Welch, K.R. & Safford, H.D. & Butz, R.J. &
O’Hara, K.L. & Stephens, S.L. (2018) From the stand scale to the landscape
scale: predicting the spatial patterns of forest regeneration after
disturbance. Ecological Applications 28(6):1626-1639.
https://doi.org/10.1002/eap.1756.

[5] Davis, K.T. & Robles, M.D. & Kemp, K.B. & Campbell, J.L (2023) Reduced
fire severity offers near-term buffer to climate-driven declines in conifer
resilience across the western United States. Proceedings of the National
Academy of Sciences of the United States of America 120(11).
https://doi.org/10.1073/pnas.2208120120.

[6] Moya, D. & González-De Vega, S. & Lozano, E. & García-Orenes, F. & Mataix-
Solera, J. & Lucas-Borja, M.E. & de las Heras, J. (2019) The burn severity and
plant recovery relationship affect the biological and chemical soil properties
of Pinus halepensis Mill. stands in the short and mid-terms after wildfire.
Journal of Environmental Management 235:250-256.
https://doi.org/10.1016/j.jenvman.2019.01.029.

[7] Scaven, V.L. & Rafferty, N.E. (2013) Physiological effects of climate
warming on flowering plants and insect pollinators and potential consequences
for their interactions. Current Zoology 59(3):418-426.
https://doi.org/10.1093/czoolo/59.3.418.

[8] Liu, Q. & Fu, B. & Chen, Z. & Chen, L. & Liu, L. & Peng, W. & Liang, Y. &
Chen, L. (2022) Evaluating effects of post-fire climate and burn severity on
the early-term regeneration of forest and shrub communities in the San Gabriel
mountains of California from Sentinel-2(MSI) images. Forests 13(7).
https://doi.org/10.3390/f13071060.

[9] Bright, B.C. & Hudak, A.T. & Kennedy, R.E. & Braaten, J.D. & Khalyani,
A.H. (2019) Examining post-fire vegetation recovery with Landsat time series
analysis in three western North American forest types. Fire Ecology 15(8).
https://doi.org/10.1186/s42408-018-0021-9.

[10] Kondylatos, S. & Prapas, I. & Ronco, M. & Papoutsis, I. & Camps-Valls, G.
& Piles, M. & Fernández-Torres, M. & Carvalhais, N. (2022) Wildfire danger
prediction and understanding with deep learning. AGU Advancing Earth and Space
Science 49(17). https://doi.org/10.1029/2022GL099368.

[11] Didan, K. (2021) MODIS/Terra vegetation indices 16-day L3 global 250m SIN
grid V061. NASA EOSDIS Land Processes Distributed Active Archive Center.
https://doi.org/10.5067/MODIS/MOD13Q1.061.

[12] Collins, C.G. & Elmendorf, S.C. & Hollister, R.D. & Henry, G.H.R & et al.
(2021) Experimental warming differentially affects vegetative and reproductive
phenology of tundra plants. Nature Communications 3442.
https://doi.org/10.1038/s41467-021-23841-2.

[13] Wan, Z. & Hook, S. & Hulley, G. (2015) MOD11A2 MODIS/Terra land surface
temperature/emissivity 8-day L3 global 1km SIN grid V006. NASA EOSDIS Land
Processes Distributed Active Archive Center.
https://doi.org/10.5067/MODIS/MOD11A2.006.

[14] Huffman, G.J. & Stocker, E.F. & Bolvin, D.T. & Neklin, E.J. & Tan, J.
(2019) GPM IMERG early precipitation L3 1 day 0.1 degree x 0.1 degree V06.
Goddard Earth Sciences Data and Information Services Center.
https://doi.org/10.5067/GPM/IMERGDE/DAY/06

[15] Shi, X. & Chen, Z & Wang, H. & Yeung, D.-Y. & Wong, W.-K. & Woo, W.-C.
(2015) Convolutional LSTM network: a machine learning approach for
precipitation nowcasting. https://doi.org/10.48550/arXiv.1506.04214

[16] LeCun, Y. & Bottou, L. & Yoshua, B. & Haffner, P. (1998) Gradient-based
learning applied to document recognition. Proceedings of the IEEE
86(11):2278-2324. https://doi.org/10.1109/5.726791.

[17] Hochreiter, S. & Schmidhuber, J. (1997) Long short-term memory. IEEE
Xplore Neural Computation 9(8):1735-1780.
https://doi.org/10.1162/neco.1997.9.8.1735.

[18] Kawano, T. & Wallbridge, N. & Plummer, C. (2020) Logistic models for
simulating the growth of plants by defining the maximum plant size as the
limit of information flow. Plant Signaling & Behavior 15(2).
https://doi.org/10.1080/15592324.2019.1709718.

[19] Kolda, T.G. & Bader, B.W. (2009) Tensor decompositions and applications.
Journal of Neuroscience 51(3):455–500. https://doi.org/10.1137/07070111X.

##  6 Appendix

###  6.1 K-Nearest Neighbor Clustering

Once each wildfire had a set of forecasted logistic growth values, we
clustered all wildfires across California to observe any possible patterns.
For this, we chose to use K-nearest neighbor (KNN) clustering enhanced with
Uniform Manifold Approximation and Projection (UMAP).

####  6.1.1 Uniform Manifold Approximation and Projection

Uniform Manifold Approximation and Projection (UMAP) is a learning technique
that can be used to visualize high-dimensional datasets in lower dimensions.
Before UMAP, one of the most widely-used techniques for data visualization was
t-distributed stochastic neighbor embedding (t-SNE), which constructs high-
dimension data in lower dimensions using Euclidean distance to determine
similarity between the dimensional maps. UMAP is functionally similar, with a
major difference being that UMAP utilizes Riemannian distance rather than
Euclidean distance. UMAP first calculates probability distributions over pairs
of objects in the high-dimension map, then assigns similar probability
distributions to a lower-dimension map of objects.

Figure 4: The results of UMAP after projecting our 5D points representing each
wildfire into a 2D plane.

To prepare our data for UMAP, we treated each wildfire as a single point in a
5D plane, where the values of five statistics collected represented their
coordinates in the plane: longitude, latitude, average predicted growth rate
$k$, average predicted carrying capacity $L$, and average precipitation during
the first 5 time steps $p$. Each set of values was also normalized to
eliminate variation in feature ranges that could introduce unwanted bias
during clustering. Then, we used UMAP to project those points in a 2D plane.
Figure 4 shows the projected results generated by UMAP. Note that the x- and
y-axes created by UMAP have no statistical significance in relation to the
variables represented in the original 5D data.

####  6.1.2 Clustering Results

The obtained UMAP points were then clustered using KNN with $3$, $5$, and $10$
clusters. The raw clustered UMAP results and the cluster labels plotted back
on a map of California are displayed in Figure 5. The size of the dots
represents the average magnitude of the $k$-values in each cluster. The
clustered labels were manually inspected, and two groups of wildfires stood
out as always being placed in the same cluster, regardless of the number of
clusters being created, indicating that the wildfires in each set were
extremely similar to each other in location and in recovery pattern.

Figure 5: The results of KNN clustering using 10 clusters, plotted on a map of
California.

The group of wildfires highlighted in red around the Counties of Los Angeles
and San Bernardino were all areas that had large positive $k$-values,
signifying that the vegetation in these areas was able to quickly recover back
to pre-wildfire levels. This cluster included the Canyon Fire 2 of 2017, which
had the highest $k$-value observed in our dataset of $0.611$. On the contrary,
the wildfires highlighted in brown around the Counties of Tuolumne and
Mariposa were all areas that had large negative $k$-values. This indicates
that the vegetation in the area continued to deteriorate , and that these
areas could need increased manual intervention to allow vegetation levels to
recover. This cluster also contained the Summit Complex Fire of 2017, which
had the lowest $k$-value observed in our study of $-1.038$.