Learning partial differential equations for biological transport models from noisy spatiotemporal data

TL;DR

This paper presents a novel neural network-based denoising approach for accurately learning biological PDE models from noisy spatiotemporal data, outperforming traditional methods.

Contribution

It introduces an ANN-based denoising method that improves PDE discovery accuracy from noisy biological data, surpassing existing denoising techniques.

Findings

ANN denoising outperforms finite differences and splines.

The method accurately recovers PDE models from noisy data.

Effective on biological transport PDEs like Fisher-KPP.

Abstract

We investigate methods for learning partial differential equation (PDE) models from spatiotemporal data under biologically realistic levels and forms of noise. Recent progress in learning PDEs from data have used sparse regression to select candidate terms from a denoised set of data, including approximated partial derivatives. We analyze the performance in utilizing previous methods to denoise data for the task of discovering the governing system of partial differential equations (PDEs). We also develop a novel methodology that uses artificial neural networks (ANNs) to denoise data and approximate partial derivatives. We test the methodology on three PDE models for biological transport, i.e., the advection-diffusion, classical Fisher-KPP, and nonlinear Fisher-KPP equations. We show that the ANN methodology outperforms previous denoising methods, including finite differences and polynomial regression splines, in the ability to accurately approximate partial derivatives and learn the correct PDE model.

Figures

Tables

Table 1. Table 1 : The relative mean-squared error (RMSE) between the noiseless data or true derivative values and our desnoised data or derivative computations using finite differences, splines, and the ANN for the diffusion-advection equation. “FD” denotes finite differences, “SP” denotes splines, and “ANN” denotes the ANN. Bold denotes the lowest errors of the three methods.

Noise	Method	$u$ RMSE	$u_t$ RMSE	$u_x$ RMSE	$u_{xx}$ RMSE
	FD	0.00e+00	5.39e-05	5.77e-04	3.69e-02
$\sigma$=0	SP	1.22e-02	1.12e+00	3.33e-02	4.86e+01
	ANN	2.86e-04	3.96e-01	8.47e-03	3.75e-01
	FD	1.02e-04	2.08e+02	4.34e-01	3.52e+01
$\sigma$=0.01	SP	1.11e-02	7.34e+00	4.93e-02	4.87e+01
	ANN	8.40e-04	7.71e-02	1.15e-02	6.93e-01
	FD	2.51e-03	2.42e+03	9.97e+00	1.01e+03
$\sigma$=0.05	SP	1.19e-02	7.47e+01	2.62e-01	5.01e+01
	ANN	5.61e-04	1.95e-01	7.71e-03	7.90e-01
	FD	1.00e-02	4.04e+03	7.59e+01	3.78e+03
$\sigma$=0.10	SP	1.04e-02	2.32e+02	1.38e+00	5.45e+01
	ANN	9.51e-04	1.23e-01	1.44e-02	7.69e-01
	FD	6.28e-02	9.04e+04	2.20e+02	3.65e+04
$\sigma$=0.25	SP	2.51e-02	6.21e+03	5.83e+00	2.42e+02
	ANN	7.29e-03	1.53e+00	4.49e-02	7.21e-01
	FD	2.41e-01	2.58e+06	1.39e+03	1.05e+05
$\sigma$=0.50	SP	3.71e-02	3.78e+04	6.68e+01	7.95e+02
	ANN	6.34e-02	3.43e+00	1.05e-01	1.44e+00

Table 2. Table 2 : Relative mean square error for Fisher-KPP equation. “FD” denotes finite differences, “SP” denotes splines, and “NN” denotes the ANN. Bold denotes the lowest errors of the three methods.

Noise	Method	$u$ RMSE	$u_t$ RMSE	$u_x$ RMSE	$u_{xx}$ RMSE
	FD	0.00e+00	3.16e-05	4.46e-04	4.67e-03
$\sigma =0$	SP	6.28e-05	2.83e-04	2.83e-03	2.01e-01
	ANN	4.86e-04	6.98e-02	1.18e-01	2.66e+00
	FD	9.82e-05	4.84e+00	6.47e+01	1.19e+03
$\sigma =0.01$	SP	6.91e-05	1.89e-01	1.56e+00	2.19e+00
	ANN	3.90e-04	3.85e-02	1.19e-02	2.55e+00
	FD	2.52e-03	9.49e+01	9.43e+02	6.78e+04
$\sigma =0.05$	SP	2.20e-04	8.85e+00	5.83e+01	2.39e+02
	ANN	4.67e-04	1.03e-02	1.63e-02	1.57e+00
	FD	9.95e-03	4.05e+02	2.81e+03	2.02e+05
$\sigma =0.10$	SP	6.80e-04	1.64e+01	6.66e+01	9.35e+02
	ANN	8.46e-04	4.26e-02	5.29e-02	1.23e+00
	FD	6.30e-02	2.92e+03	3.29e+04	3.92e+05
$\sigma$=0.25	SP	4.14e-03	2.26e+02	5.65e+02	1.44e+04
	ANN	6.41e-03	6.94e-02	1.04e-01	5.90e+00
	FD	2.38e-01	9.22e+03	1.01e+05	1.04e+06
$\sigma =0.50$	SP	1.52e-02	5.44e+02	1.32e+03	3.22e+04
	ANN	6.60e-02	5.48e-01	8.98e-01	3.52e+01

Table 3. Table 3 : Relative mean square error for nonlinear Fisher-KPP equation. “FD” denotes finite differences, “SP” denotes splines, and “ANN” denotes the ANN. Bold denotes the lowest errors of the three methods.

Noise	Method	$u$ RMSE	$u_t$ RMSE	$u_x$ RMSE	$u_{xx}$ RMSE
	FD	9.71e-36	2.20e-05	1.53e-04	1.34e-01
$\sigma =0$	SP	1.53e-05	4.19e-05	1.09e-03	5.71e+00
	ANN	8.73e-04	1.41e+01	6.16e+00	1.40e+02
	FD	1.02e-04	2.21e+02	9.16e+03	8.26e+02
$\sigma =0.01$	SP	2.01e-05	1.02e+01	6.20e+02	3.32e+00
	ANN	6.87e-04	5.55e+00	6.58e+01	1.90e+02
	FD	2.43e-03	5.45e+03	2.65e+05	3.35e+04
$\sigma =0.05$	SP	1.57e-04	2.74e+02	6.84e+03	6.69e+01
	ANN	1.08e-03	8.26e+00	1.68e+00	1.70e+02
	FD	1.01e-02	2.30e+04	1.16e+06	5.69e+04
$\sigma =0.10$	SP	6.20e-04	1.32e+03	1.73e+04	7.23e+01
	ANN	1.84e-03	1.93e+01	1.10e+01	2.04e+02
	FD	6.25e-02	1.47e+05	4.95e+06	9.86e+05
$\sigma =0.25$	SP	4.00e-03	5.69e+03	1.24e+05	1.54e+03
	ANN	6.34e-03	2.58e+01	2.34e+01	2.59e+02
	FD	2.43e-01	5.73e+05	1.92e+07	2.46e+06
$\sigma =0.50$	SP	1.38e-02	2.87e+04	4.16e+05	1.13e+04
	ANN	6.89e-02	6.49e+01	1.97e+02	4.88e+02

Table 4. Table 4 : Inferred parameters for the nonlinear Fisher-KPP Equation data when performing an inverse problem on Equation ( 18 ).

	True Equation
	$u_t=.02u{u}_{xx}+.02u_x^2+10u-10u^2$
$\sigma$	Revised Equation
0	$u_t=2.3\times {10}^{-12}{u}_{xx}+.017u{u}_{xx}+.021u_x^2+10.0u-10.0u^2-1.60\times {10}^{-8}$
0.01	$u_t=1.4\times {10}^{-5}{u}_{xx}+.021u{u}_{xx}+.019u_x^2+10.0u-10.0u^2-1.38\times {10}^{-8}$
0.05	$u_t=2.5\times {10}^{-4}{u}_{xx}+.0034u{u}_{xx}+.029u_x^2+9.9u-10.0u^2+-2.86\times {10}^{-8}$
0.10	$u_t=6.1\times {10}^{-4}{u}_{xx}+8.42\times {10}^{-4}u{u}_{xx}+.023u_x^2+9.38u-9.52u^2-2.44\times {10}^{-4}$
0.25	$u_t=7.2\times {10}^{-4}{u}_{xx}+.015u{u}_{xx}+.024u_x^2+9.8u-9.14u^2+3.53\times {10}^{-7}$
0.50	$u_t=2.1\times {10}^{-3}{u}_{xx}-3.72\times {10}^{-3}u{u}_{xx}+.032u_x^2+9.7u-7.5u^2+1.38\times {10}^{-6}$

Equations

$$u_t(x,t)$$

$$u(x,t_0)$$

$$U_{i,j}$$

$$w_{i,j}=\left({\beta }_1{u}_1^{\gamma }(x_i,t_j),\dots ,{\beta }_n{u}_n^{\gamma }(x_i,t_j)\right){}^{\text{T}}$$

$$u_t=-c{u}_x+D{u}_{xx},D,c\in \mathbb{R}$$

$$u_t=D{u}_{xx}+ru-r{u}^2,D,r\in \mathbb{R}$$

$$u_t=Du{u}_{xx}+D{u}_x^2+ru-r{u}^2,D,r\in \mathbb{R}$$

$$h_{\theta}(\vec{x})=a_{2}\bigg{(}W_{2}\Big{(}a_{1}(W_{1}\vec{x}+b_{1})\Big{)}+b_{2}\bigg{)}$$

$$\mathcal{J}(\theta )=\frac{1}{MN}\sum _{i=1,j=1}^{M,N}{\left(\frac{h_{\theta }(x_i,t_j)-u_{i,j}}{|h_{\theta }(x_i,t_j){|}^{\gamma }}\right)}^2.$$

$$\mathcal{L}(\theta )=\frac{1}{MN}\sum _{i=1,j=1}^{M,N}\left[{\left(\frac{h_{\theta }(x_i,t_j)-u_{i,j}}{|h_{\theta }(x_i,t_j){|}^{\gamma }}\right)}^2+\lambda {\left(W_1{\vec{x}}_{i,j}+b_1\right)}^2\right]+\frac{1}{MN}\sum _{h_{\theta }\notin [0,1]}{h}_{\theta }^2$$

$$\Theta =\text{setcounter}MaxMatrixCols13\left[\begin{array}{ccccccccccccc}1& u& \cdots & u^p& u_x& \cdots & u^p\odot u_x& u_{xx}& \cdots & u^p\odot u_{xx}& u_x^2& u_x\odot u_{xx}& u_{xx}^2\end{array}\right]$$

$$u_t\approx \Theta \xi ,$$

$$\hat{\xi }=\arg \underset{\xi \in {\mathbb{R}}^d}{\min }\frac{1}{MN}\Vert u_t-\Theta \xi {\Vert }_2^2,\text{subject to }\Vert \xi {\Vert }_0\le k$$

$$u_t^{\text{train}}={\tilde{\Theta }}_i^{\text{train}}{\xi }_i.$$

$$\text{TPR}(\xi )=\frac{\text{TP}}{\text{TP + FN + FP}},$$

$$u_t=a+b{u}_x+c{u}_{xx}+d{u}^2+eu+f{u}^2{u}_x+g{u}^2{u}_{xx},a,...,g\in \mathbb{R}.$$

$$u_t=a{u}_x+b{u}_{xx}+c{u}^2{u}_x,a,b,c\in \mathbb{R}$$

$$u_t=a{u}_x+b{u}_{xx}+c{u}^2{u}_x+d{u}^2{u}_{xx}+e{u}_x^2,a,...,e\in \mathbb{R},$$

$$u_t=a{u}_{xx}+bu{u}_{xx}+c{u}_x^2+du+e{u}^2+f,a,\dots ,f\in \mathbb{R}.$$

Full text

Keywords: Interpolation, numerical Differentiation, equation learning, sparse regression, partial differential equations, parameter estimation

1 Introduction

Recent research has investigated methods for discovering systems of differential equations that describe the underlying dynamics of spatiotemporal data. There are key advantages to learning and then using mathematical models for prediction instead of using a purely machine learning based method, e.g, neural networks. First, if the learned mathematical model is an accurate description of the processes governing the observed data, it has the ability to generalize from the set of training data to data outside of the training domain. Second, the learned mathematical model is interpretable, making it informative for scientists to hypothesize the underlying physical or biological laws governing the observed data. Examples of recent methods for inferring the underlying governing equations include the Sparse Identification of Nonlinear Dynamics (SINDy) algorithm [10] and the Equation Learner (EQL) neural network [23, 30], both of which are used for discovering systems of ordinary differential equations (ODEs), and the PDE Functional Identification of Nonlinear Dynamics (PDE-FIND) algorithm [28], which is used to identify PDE systems. Boninsegna et. al [8] recently extended the SINDy algorithm to recover stochastic dynamical systems. Model selection criteria (such as Akaike Information Criteria and Bayesian Information Criteria) have been combined with the SINDy algorithm to increase robustness to errors, although incorrect models were still selected at noise levels we consider here [22]. The discovery methods mentioned above assume that the measured data arise from a parameterized $n$-dimensional dynamical system of the form

[TABLE]

with parameter vector $\theta\in\mathbb{R}^{k}$ and appropriate boundary conditions. It is then assumed that the true dynamical system, $F(\cdot)$, is only comprised of a few terms (e.g., for the diffusion-advection equation, $F=Du_{xx}-cu_{x},D,c\in\mathbb{R}$). The goal of these methods is to correctly specify the small number of correct terms from a large library with many potential candidate terms (such as $u$, $u_{x}$, $u_{xx}$, $u^{2}$, $uu_{x}$, etc.) with the aid of sparse regression [26]. A practical challenge arises because one typically does not have access to the noiseless values of $u(x,t)$ or its partial derivatives. Instead, one must approximate these values from noisy experimental data. Here, our goal is to investigate the performance of existing denoising methods that are used in conjunction with PDE-FIND and to present a novel denoising methodology relying on artificial neural networks (ANNs).

Several methods have been used for denoising data to approximate $u(x,t)$ and it partial derivatives ($u_{t}$, $u_{x}$, $u_{xx}$, etc.). The most prevalent methods that have been proposed are finite difference approximations or the use of cubic splines for interpolation, followed by partial differentiation of the fitted splines. However, both of these methods have been found to be prone to inaccuracies in the presence of noise [28]. Recent work has considered recovery of dynamical systems with high amounts of noise added to the time derivative measurement ($u_{t}$) by transforming the data into a spectral domain [31]. Zhang et al. [35] recently proposed using sparse Bayesian regression, which allows for error bars for each candidate term in the discovered equation. However, although their method was robust, the noise in this study was also added only to the time derivative term ($u_{t}$) instead of the observed data ($u(x,t)$). Importantly, it has been noted that introducing noise to the observed data itself ($u(x,t)$) hinders the recovery of the correct PDE, thus, developing a method of denoising data for $u(x,t)$ has been identified as a current challenge for learning PDEs [31]. To the best of our knowledge, the use of finite differences or utilizing splines are the two methods that, in practice, yield the most accurate approximations for the library terms involved in PDE learning. The primary challenge involved with using these methods for numerical differentiation, which we further test in this work, is that they are sensitive to noise levels and can amplify noise as the order of the derivative increases. This challenge inhibits learning PDEs for practical biological applications where data may have large noise levels due to many sources of error, including the data collection process, imprecise measurement tools, and the inherent stochastic nature of biological processes [5, 12]. For example, for ecological measurements of population abundance, typical data sets can have noise levels on the order of a coefficient of variation equal to 0.2 [13, 25]. Notably, adding this biologically relevant level of noise to the observation $u(x,t)$ has not been considered in previous PDE learning work.

An additional, yet realistic, complication that has not been considered is the presence of non-constant error noise in the spatiotemporal data used for PDE learning. For example, proportional error noise can occur when the variance of the data is proportional to the size of the measurement, e.g., population size or density [6]. Non-constant error noise may also occur when the observed processes occur on different time scales [17]. To account for non-constant error noise in the scenario that one has a mathematical model for the biological process generating the data, e.g., $u(x,t)$, the non-constant error noise can be accounted for with a statistical model used in conjunction with the mathematical model [2]. For example, for a set of observed data at space points $x_{i}$, $i=1,\dots,M$ and time points $t_{j}$, $j=1,\dots,N$, a general statistical model is given by

[TABLE]

where the noiseless observations are corrupted by noise modeled by the random variable $w_{i,j}\odot\mathcal{E}_{i,j}$ in which $\odot$ represents element-wise multiplication. Finite difference methods assume $w_{i,j}\odot\mathcal{E}_{i,j}=0$ while regression methods using splines assume the variance of $w_{i,j}\odot\mathcal{E}_{i,j}$ is constant. More generally, the error term $\mathcal{E}_{i,j}$ may instead be generated by a probability distribution that is weighted by

[TABLE]

for $\gamma\geq 0$ and $\beta_{1},\dots,\beta_{n}\in\mathbb{R}$. Constant error noise is modeled by assuming $\gamma=0$ and $\beta_{1},\dots,\beta_{n}=1$. Proportional error noise is modeled by assuming $\gamma>0$, $\beta_{1},\dots,\beta_{n}\neq 0$ [7]. We note that finite difference and spline approximations of noisy data neglect whether the noise process has non-constant variance. We hypothesize that this leads finite difference and spline approximations to yield poor estimates of the noiseless data $u(x,t)$ and its partial derivatives when the data contain proportional error noise. In this work, we investigate this hypothesis and develop a methodology using ANNs as a model for $u(x,t)$ in conjunction with an appropriate statistical error model that accounts for the presence of proportional error when denoising spatiotemporal data.

The denoising methods present in this work focus on spatiotemporal data for learning PDEs, however the methods we describe can be readily applied to learning ODEs. We choose to focus our study on a specific set of diffusive PDE models, which have provided a wealth of insight into many biological transport phenomena, including ecological migration and invasion [14], neuronal transport [18], cancer progression [3, 34, 27, 29], and wound healing [21, 24]. Here, we demonstrate how an ANN can be used with a non-constant error statistical model to accurately approximate $u(x,t)$ from noisy proportional error data. It has long been known that ANNs are universal function approximators [15], meaning ANNs have the capacity to approximate continuous functions arbitrarily well. Unlike local approximations such as finite differences and splines, an ANN can be fit to an entire set of spatiotemporal data, affording a global context that may help to decrease overfitting. We compare the accuracy of the ANN-based method to finite differences and splines in computing estimates for $u(x,t)$ and its partial derivatives. We then investigate whether the PDE-FIND algorithm can more accurately infer the governing PDE equations from data when its library of terms is constructed using an ANN-based method.

2 Methods

The process of learning a system of equations from noisy data can be divided into two main components: (1) the data denoising and library construction component, in which the underlying dynamical system $u(x,t)$ and its partial derivatives are approximated from the noisy realizations of \big{\{}U_{i,j}\big{\}}_{i=1,j=1}^{M,N} in (2), and (2) the equation learning component, in which, given approximations for $u$, $u_{t}$, $u_{x}$, $u_{xx}$, etc., one employs an algorithm that can effectively uncover the mechanistic form of $F$ in (1a) (Figure 1). Below, we describe the mathematical models used for data generation, the method of constructing a library from noisy data, and equation learning. All of the denoising methods were implemented in Python 2.7 using the Scipy package for polynomial splines and the Keras machine learning library for ANNs. All code and videos are available at https://github.com/biomathlab/PDElearning/.

2.1 Data Generation

We consider three diffusive PDE models for biological transport in this work, each of which has been used previously to interpret biological data [16, 21, 32, 33]. These models include the diffusion-advection equation

[TABLE]

the classical Fisher-Kolmogorov-Petrovsky-Piskunov (Fisher-KPP) Equation

[TABLE]

and the nonlinear Fisher-KPP Equation

[TABLE]

where $D$ is the diffusion coefficient, $r$ is the intrinsic population growth rate, and $c$ is the advection rate.

Assume $u(x,t)$ denotes the solution to one of the above mathematical models. We generate noisy data by using Equation (2) with $w_{i,j}=\sigma u(x_{i},t_{j})$ (i.e., $\beta=\sigma\text{ and }\gamma=1$) in which all $\mathcal{E}_{i,j}$ terms are simulated as i.i.d. normal random variables with mean zero and variance one. We generate six data sets for each mathematical model, setting $\sigma=$ 0, 0.01, 0.05, 0.10, 0.25, and 0.50. For Equation (4), we use its analytical solution to compute $u(x,t)$. For Equations (5) and (6) we use finite difference computations to numerically approximate $u(x,t)$. The numerical step sizes in these computations were chosen small enough to not introduce significant noise into the solution. We use $M=101$ spatial locations and $N=300$ time points to generate data for the diffusion-advection equation and $M=199$ spatial locations and $N=99$ time points for the Fisher-KPP and nonlinear Fisher-KPP Equations. As a preprocessing step, each data set is scaled to $[0,1]$ in order to consistently measure errors across the various data sets and noise levels.

2.2 Data denoising and library construction

2.2.1 Finite differences and spline approximations

For finite differences, we use central difference formulas on interior points and forward differences at the boundaries to obtain first order derivative approximations. For higher order derivatives (e.g. $u_{xx}$), first order finite difference rules are repeated on the corresponding previous order derivative approximations.

Finite difference approximations can be obtained accurately, efficiently, and directly from noiseless data, however, their accuracy quickly deteriorates in the presence of observation error. Following [28], we also employ polynomial spline regression. For a given data point $u_{i,j}$ in the set of observations, we fit a cubic bi-spline on a small two-dimensional neighborhood of size $11\times 11$ centered at $u_{i,j}$. Denoised function and derivative approximations are then obtained using evaluations and the analytic derivatives of the fitted polynomial at the center point $u_{i,j}$. Note that since we only approximate derivatives up to second order, higher order splines are not considered. For points close to the boundaries, we evaluate along the spline approximation that is nearest to the boundary. Each polynomial spline is optimized using ordinary least squares. We found that cubic bi-spline approximations were generally more robust than the one-dimensional cubic splines used in [28] (Supplementary Figures S1-S3).

2.2.2 Artificial neural network approximations

An artificial neural network (ANN), denoted $h_{\theta}(\vec{x})$, was used to approximate $u(x,t)$, with one hidden layer of the form

[TABLE]

where $\vec{x}=[x\,\,\,t]{{}^{\text{T}}}$ and $a_{i}(\cdot)$ represents the continuous nonlinear activation function for the i${}^{\text{th}}$ layer. The matrices $W_{i}$ and vectors $b_{i}$ (typically called weights and biases) comprise the total set of trainable parameters $\theta=\{W_{1},b_{1},W_{2},b_{2}\}$ of the network. We note that the use of one hidden layer in a neural network is a sufficient condition to make it a universal function approximator under the assumption that the activation function is bounded and non-constant [15]. This result extends to ANNs with multiple hidden layers, however, we found that while training multi-layer ANNs resulted in faster convergence, the derivative approximations were worse. The task is therefore to find the optimal parameters $\theta^{*}$ such that $h_{\theta^{*}}(x,t)\approx u(x,t)$. The fitted surface function $h_{\theta^{*}}(x,t)$ and the computation of analytic derivatives of this function are used to approximate $u(x,t)$ in (2) and its partial derivatives for library construction in the PDE learning task.

We formulate a regression problem using the ANN defined in (7) in order to find $\theta^{*}$ that minimizes the generalized least squares cost function

[TABLE]

Note that this error formulation accounts for the statistical error model in (2) and also reduces to ordinary least squares when $\gamma=0$. Also note that function values $h_{\theta}(x_{i},t_{j})$ less than 1e-4 in absolute value are set equal to one in the denominator $|h_{\theta}(x_{i},t_{j})|^{\gamma}$ during network training and evaluation. The network parameters $\theta$ are optimized using the first-order gradient-based “Adam” optimizer [19] with default parameters and a batch size of 10.

We found that the choice of activation function in the ANN, $a_{i}(\cdot)$, plays an important role in the accuracy of the partial derivative approximations. Typical activations like sigmoid and hyperbolic tangent yield oscillations in higher order derivative terms (e.g., see Supplementary movie S1.). To mitigate this, we chose to use the “softplus” activation function which takes the form $\log(1+e^{z})$. This function has many desirable properties for approximating $u(x,t)$ (e.g. smoothness and infinitely many derivatives), however, it is unbounded, which violates an assumption of ANNs as universal approximators. We address the unboundedness of the softplus function by including an $\ell_{2}$-regularization penalty on the inputs $W_{1}\vec{x}+b_{1}$ in (7). Further, an additional squared error term is included in the loss function to penalize function values outside $[0,1]$. Without this term, the function values can blow up during training since $\mathcal{J}(\theta)\rightarrow 1$ as $h_{\theta}(\vec{x})\rightarrow+\infty$ when $\gamma=1$. Thus, the complete loss function used for training the ANN is

[TABLE]

where the first term corresponds to the generalized least squares cost function (8), the second term corresponds to the $\ell_{2}$-regularization penalty with regularization strength $\lambda$, and the third term corresponds to the additional squared error term to penalize function values outside $[0,1]$. Generally, the regularization strength should be chosen separately for each data set using cross-validation, but we found that no regularization (i.e. $\lambda=0$) was needed for the 18 data sets considered in this paper.

We used 1,000 neurons in the hidden layer of the ANN. This choice was large enough to have maximal capacity to fit the data, while still allowing the optimization of $\theta$ to be computationally feasible on a desktop computer (3.4 GHz Intel Core i5 processor, 8gb RAM). In order to prevent overfitting, the data were randomly split into training and validation sets (90%/10%) when training the ANN. The optimal network parameters were chosen to be those that minimized the error (8) on the validation set. Training ended if the validation error did not decrease for 50 consecutive epochs.

Representative examples of results from the bi-spline and ANN methods are shown in Supplementary movies S2 and S3, respectively. All movies for all methods and noise levels considered in this work can be found at https://github.com/biomathlab/PDElearning/animations/.

2.3 Equation Learning

We use the PDE-FIND algorithm [28] to discover the form of $F(\cdot)$ in Equation (1a) using computations of $u,u_{x},u_{xx},$ and $u_{t}$ from the ANN, spline, and finite difference methods. Prior to implementing PDE-FIND, the numerical approximations are scaled from $[0,1]$ back into their original scales. We discuss the PDE-FIND implementation in Section 2.3.1 and an additional pruning method in Section 2.3.2 that is used to remove extra terms from the final learned equation. We further discuss how we analyze our results in Section 2.3.3.

2.3.1 PDE-FIND Implementation

Once $u(x,t)$ and its partial derivatives have been computed, a large library of potential PDE terms is formed column-wise in the matrix, $\Theta$, given by

[TABLE]

where each column of $\Theta$ is some vectorization of the written term. The ANN has difficulty capturing the early dynamics of the diffusion-advection equation, so we skip the first 20 timepoints from the denoised data when building $\Theta$ for all data sets and denoising strategies. To reduce the computational time, only every fifth remaining timepoint is included in $\Theta$. Hence, while the data sets for the diffusion-advection equation begin with $N=300$ timepoints, only (300-20)/5=56 timepoints are used in constructing $\Theta$. We set $p=2$ resulting in $d=12$ columns in $\Theta$. Each column of $\Theta$ thus represents a candidate term comprising $F$, so we assume

[TABLE]

where $\xi$ is a vector whose nonzero entries correspond to the true terms of $F$. The vector $\xi$ is estimated using methods from sparse regression [26]. Sequential threshold ridge regression was found to be a suitable method for estimating $\xi$ for PDE-FIND in a previous study [28]. However, we found that the Greedy algorithm to performed well for the data and models we considered in this work. The Greedy algorithm computes

[TABLE]

for some sparsity parameter, $k$ [36].

The value of $k$ for a given data set is treated as a hyperparameter that is found by splitting the library data into separate training and validation sets, and then optimizing over the validation set. In this training-validation procedure, we randomly divide our data points for $u_{t}$ into 5-by-5 tiles of adjacent spatiotemporal points and then randomly assign 50% of these tiles to a training data set, $u_{t}^{\text{train}}$, and the remaining 50% to a validation set, $u_{t}^{\text{validate}}$. We split the corresponding rows of $\Theta$ into $\Theta^{\text{train}}$ and $\Theta^{\text{validate}}$. We perform our hyperparameter search over 51 values for $k$ between 0 and $10^{3}$. For each value of $k$, we estimate $\hat{\xi}$ from the training set. For each estimate, we then compute its mean-squared error (MSE) over the validation set. We choose the hyperparameter corresponding to the $\hat{\xi}$ estimate with the smallest MSE on the validation data. The equation that results from sparse regression with this hyperparameter is our final equation from the PDE-FIND algorithm. We refer to the validation MSE from the final equation “$\text{val}_{0}$” in the remaining text.

2.3.2 Pruning Method

We chose a 50-50 training and validation split for the data to avoid overfitting to the training data with a large validation set. Even so, we will demonstrate in Section 3.2 below that PDE-FIND is able to learn small but systematic biases from the ANN’s fit to $u$ and its derivatives by incorporating extra terms into the final equations. Pruning methods have previously been developed that remove extra terms that do not significantly increase an algorithm’s performance, see for example [1, 20]. Accordingly, we implement the following pruning method after the PDE-FIND implementation described in Section 2.3.1 for all methods in order to delete the extra terms from the final equation.

The pruning procedure starts with the reduced library of candidate terms (i.e., the columns of $\Theta$) for the right hand side of Equation (1a) with nonzero entries of $\hat{\xi}$ that resulted from our training-validation procedure. We then perform a sensitivity test for the remaining terms as follows. Suppose $\tilde{d}$ terms remain in $\tilde{\Theta}$, and let $\tilde{\Theta}_{i},i=1,\ldots,\tilde{d}$ denote the further-reduced library where the $i$th column of $\tilde{\Theta}$ has been removed. For each value of $i$, we find the least squares solution (without regularization) on the training data to the equation

[TABLE]

We then use our $\hat{\xi}_{i}$ estimate and compute the MSE over the validation data when the $i^{\text{th}}$ term has been removed and call this computation $\text{val}_{i}$. We then remove any candidate terms for our library that result in $\text{val}_{i}/\text{val}_{0}<1+\alpha$ for some $\alpha>0$. After this pruning step, we perform one final round of training without regularization over the fully reduced library to find the final form of our underlying equation.

It is important to note that choice of the $\alpha$ pruning threshold value warrants careful decision. If this value is chosen too high, then too few terms will be selected and the learned equation will be incomplete. If the chosen value is too small, then the final equation will admit extra terms arising from the systematic errors in derivative estimation. We will demonstrate below that the choice of arbitrary choice of $\alpha=0.25$ below provides promising results for the diffusion-advection and Fisher-KPP Equations, while $\alpha=0.05$ is suitable for the nonlinear Fisher-KPP Equation.

2.3.3 Accuracy Metrics

To quantitatively assess the accuracy in recovering the correct PDE that generated the data, i.e., using the combined PDE-FIND with pruning methodology described above, we introduce the the true positive ratio (TPR) for a given vector $\xi$ as:

[TABLE]

where TP (“True Positive”) denotes the number of correctly-specified nonzero coefficients in $\xi$, FN (“False Negative”) denotes the number of coefficients in $\xi$ that are incorrectly specified as zero, and FP (“False Positive”) denotes the number of coefficients in $\xi$ that are incorrectly specified as nonzero. Recall that the nonzero entries of $\xi$ correspond to the relevant terms in an equation (i.e., for a library of $\Theta=[1\ u\ u_{x}\ u_{xx}],\xi=[0\ 1\ 2\ 0]^{T}$ corresponds to the equation $u_{t}=u+2u_{x}$). For example, when trying to learn Equation (4), an equation of the form $u_{t}=u_{xx}+uu_{x}$ would have TP = 1 (the nonzero coefficient for $u_{xx}$ is correct), FN = 1 (the missing $u_{x}$ term is incorrect), FP = 1 (the nonzero $uu_{x}$ term is incorrect), resulting in a final score of TPR = 1/3. Note that the TPR value is similar in nature to the Jaccard index: larger TPR values suggest that the true equation form has been better approximated, and TPR = 1 signifies that the correct equation form has been recovered.

We note that the learned equation from the PDE-FIND with pruning method was often found to be sensitive to the random split of $u_{t}$ and $\Theta$ into training and validation data. Therefore, we performed PDE learning for 1,000 different random training-validation data splits of $u_{t}$ and $\Theta$ for each data set and for each computational method (finite differences, splines, and ANN). We then consider the distribution of TPR($\hat{\xi}$) scores to assess the overall performance of the methodology. We declare the most commonly-learned equation among the 1,000 data splits as the final learned equation for each data set and computational method.

3 Results

In this section, we detail our results using the ANN to denoise data for $u(x,t)$ and compute partial derivatives. In addition, we test the accuracy of using the ANN method in conjunction with PDE-FIND to learn PDEs. Analogous results are presented for finite differences and splines. We begin by demonstrating the accuracy of the partial derivative calculations in Section 3.1, we explain why PDE-FIND finds small systematic bias terms in Section 3.2, and then we detail the accuracy in learning of the diffusion-advection, Fisher-KPP, and nonlinear Fisher-KPP equations in Sections 3.3-3.5.

3.1 Derivative Calculations

We found that the finite difference method most accurately approximates $u$ and its derivatives for $\sigma=0$ (Table 1). This result is not surprising, as finite difference computations assume that there is no error in the data. For all other values of $\sigma$, we observe that the ANN produces the most accurate derivative calculations, though the splines appear best at inferring $u(x,t)$ when the data are noisy. It is important to note that the ANN’s calculations are often several magnitudes of order more accurate than the spline and finite difference approximations, and this disparity between the computations appears to increase with $\sigma$ (Table 1). For example, at $\sigma=0.01$, the ANN’s relative mean squared error (RMSE) for $u_{t}$ is four orders of magnitude smaller than the RMSE for finite differences and two orders of magnitude smaller than the RMSE for splines. At $\sigma=0.50$, the ANN’s RMSE for $u_{t}$ has become six orders of magnitude smaller than the RMSE for finite differences and four orders of magnitude smaller than the RMSE for splines. The other derivative computations show similar results.

Similarly, we found that the ANN was most accurate for noisy data from the Fisher-KPP Equation (Table 2) and the nonlinear Fisher-KPP Equation (Table 3). Recall that we do not have analytical solutions to these equations, so we used finite difference computations on the noiseless data $(\sigma=0)$ as an estimate for the analytical derivative values for the Fisher-KPP and nonlinear Fisher-KPP Equations. Again, in both cases we observe that the finite difference calculations perform best in computing the RMSE for $\sigma=0$, but on average, the ANN provides the best calculations for the derivatives for larger values of $\sigma$. The spline method is consistently the most accurate at inferring $u$ from the data. The disparity between the RMSE calculations for the ANN as compared to the splines or finite differences again appears to increase with $\sigma$ for these two equations.

3.2 PDE-FIND without pruning learns the wrong equation

We found that, in general, the PDE-FIND method learns the wrong equation, even when no noise is added to the data (Supplementary section S2, Figures S4-S6). Each denoising method resulted in accurate estimates for $u(x,t)$ and its partial derivatives in this case, however. For example, the residuals between the ANN model and the analytical values for $u,u_{t},u_{x},\text{ and }u_{xx}$ were small when $\sigma=0$ (Figure 2). We observed that, while small, the ANN residuals include systematic biases comprised of regions of over- and under-prediction. For example, all points near $(x,t)=(0.6,0.4)$ for the ANN’s calculation for $u_{x}$ appear to over-predict the true value for $u_{x}$ in this region (Figure 2). This contradicts the assumption of independence in $\{\epsilon_{i,j}\}_{i=1,j=1}^{M,N}$ for the statistical model in Equation (2). As we will now demonstrate, these small, systematic error terms from the ANN cause PDE-FIND to learn the incorrect equation.

We illustrate here that PDE-FIND learns the incorrect equation when training data for $u_{t}$ is comprised of $u_{t}$ at all spatial points for the first half of the given time points and the validation data is comprised of all spatial points for the second half of all time points. Recall that in our actual implementation discussed below, we randomly split the training and validation data in $5\times 5$ bins of adjacent spatiotemportal points. Using denoised values for $u(x,t)$ and its partial derivatives from the ANN in the case where $\sigma=0$ in the data, our training-validation procedure without pruning learns an equation of the form

[TABLE]

Similarly, the learned equations using finite difference and spline computations are

[TABLE]

and

[TABLE]

respectively.

Each of these equations are incorrect and have extra terms on the right hand side of the learned PDE for the diffusion-advection equation. In Figure 3, we depict illustrative portions of the training and validation sets comparing the analytical values of $u_{t}$ against the computed values of $u_{t}$ and PDE-FIND’s selected equation using ANN approximations. We found that PDE-FIND selects Equation (15) in place of the true diffusion-advection equation because it recovers the ANN’s incorrect computations of $u_{t}$ in both the training and validation data. In doing so, PDE-FIND fits the erroneous $u_{t}$ computations from the ANN approximation by including extra terms in the learned PDE.

3.3 PDE-FIND with pruning for the diffusion-advection equation

We tested whether implementing an additional pruning step with PDE-FIND could remove the extra terms resulting from the biases discussed in Section 3.2. For the diffusion-advection equation, we found that for all values of $\sigma$ except $\sigma=0.01$, PDE-FIND with pruning achieves the highest median TPR when using ANN approximations (Figure 4). The ANN’s median value is TPR = 1 (meaning that over half of the simulations yielded the correct equation form) for $\sigma=0,0.05,0.10,\text{ and }0.25$. The ANN resulted in a median TPR = 0.667 at $\sigma=0.50$. In contrast, the spline method only achieved of median TPR = 1 at the lower noise levels $\sigma=0.01\text{ and }0.05$. For $\sigma\geq 0.10$, the medians for the spline method were all TPR$\leq$0.5. The finite difference method resulted in a median TPR = 1 at $\sigma=0.01$, but the median TPR = 0 for larger values of $\sigma$.

Supplementary Table S1 shows the most commonly learned PDEs for each denoising method at each noise level. We found that the ANN method, used in conjunction with PDE-FIND with pruning, resulted in the correct PDE for $\sigma=0,.05,0.10,0.25$. The ANN specifies the incorrect equation for $\sigma=0.01$ and $\sigma=0.50$. However in both of these cases, the extra terms have small parameter values (e.g. 0.001) that a scientist with an understanding of the system under consideration may manually neglect. On the other hand, PDE-FIND cannot discover the correct equation with finite difference or spline computations for $\sigma\geq 0.10$. These results suggest that the ANN method enables PDE-FIND with pruning to learn the diffusion-advection equation accurately at biologically realistic noise levels, e.g, $\sigma=.05,0.10,0.25$.

3.4 PDE-FIND with pruning for the Fisher-KPP Equation

We tested the PDE-FIND with pruning method in conjuction with several denoising strategies using data from the Fisher-KPP Equation. We found that the ANN method had a median TPR = 1.0 (meaning that the correct equation is specified for at least half of the training-validation data splits) for $\sigma=0,0.01,0.05,\text{ and }0.10$ (Figure 5). In contrast, the finite difference calculations only had a median TPR = 1 at $\sigma=0$, and the spline method only had median TPR = 1 at $\sigma=0,0.01$. The accuracy in using PDE-FIND with the spline and finite difference methods quickly deteriorates for high noise levels. The finite difference method resulted in a median TPR = 0 for $\sigma\geq 0.05$, and the spline method resulted in a median TPR = 0.667 at $\sigma=0.05,0.10$ and TPR = 0 at $\sigma=0.25,0.50$. The ANN had a median TPR = 0.6 and 0.5 for $\sigma=0.25\text{ and }0.5,$ respectively.

Supplementary Table S2 shows the most commonly chosen PDEs resulting from the PDE-FIND with pruning method. We found that PDE-FIND with pruning is able to discover the correct equation form for $\sigma=0,0.01,0.05,\text{ and }0.10$ when using the ANN approximations. While PDE-FIND is unable to specify the correct equations with ANN data for $\sigma=0.25$ and 0.50, all of the terms in the Fisher-KPP equation were included in the learned PDEs. In contrast, using the spline method for denoising resulted in only learning the correct PDE for $\sigma\leq 0.01$. For larger values of $\sigma$, the spline method resulted in large errors in the derivative approximations and did not yield any terms on the right hand side of the learned PDE. Similarly, the finite difference method resulted in only learning the true equation form for $\sigma=0$. These results suggest that only the ANN method enables PDE-FIND with pruning to learn the Fisher-KPP equation for reasonably high noise levels of $\sigma=0.05,0.10$.

3.5 PDE-FIND with pruning for the nonlinear Fisher-KPP Equation

We found that the PDE-FIND with pruning method struggles with all three denoising strategies to recover the correct PDE from data that has been generated by the nonlinear Fisher-KPP Equation. PDE-FIND could not achieve a median TPR = 1 for any of these methods, meaning that the correct equation was never specified for over half of the training-validation data splits (Figure 6). All methods have median TPR = 0.8 at $\sigma=0$. When using ANN approximations, the PDE-FIND with pruning method has median TPR = 0.8 at $\sigma=0.01\text{ and }0.05$, TPR = 0.6 at $\sigma=0.10\text{ and }0.25$, and TPR = 0.5 at $\sigma=0.50$. When using spline computations, PDE-FIND with pruning has median TPR = 0.8 at $\sigma=0.01$, TPR = 0.5 at $\sigma=0.05,0.10,\text{ and }0.25$, and TPR = 0 at $\sigma=0.50$. When using finite difference computations, PDE-FIND with pruning has median = 0.50 at $\sigma=0.01\text{ and }0.05$ and TPR = 0 for $\sigma\geq 0.10$.

While all of the denoising strategies lead to incorrect equations, the ANN strategy recovers the most relevant terms in its final equations. Supplementary Table S3 shows the most-commonly chosen PDEs resulting from the PDE-FIND with pruning method. We found that all three methods predict the true equation with an extra Fickian diffusion term, $u_{xx}$, for $\sigma=0$. When using the ANN approximations, the PDE-FIND with pruning algorithm recovers the true equation with added Fickian diffusion at $\sigma=0.05$, and it recovers three of the correct terms but excludes $uu_{xx}$ at $\sigma=0.01$. For larger values of $\sigma$ with the ANN approximations, PDE-FIND with pruning recovers three correct terms, excludes the $uu_{xx}$ term, and includes an extra Fickian diffusion term (as well as an additional constant term at $\sigma=0.50$). When using spline computations, the PDE-FIND with pruning algorithm recovers three correct terms, excludes the $uu_{xx}$ term, and adds an extra Fickian diffusion term at $\sigma=0.01$. For $\sigma=0.05-0.25$, the final equation recovers two correct terms but excludes the $uu_{xx}$ and $u_{x}^{2}$ terms. At $\sigma=0.50,$ all terms are deleted when using spline computations. When using finite difference computations, the PDE-FIND with pruning algorithm correctly recovers two terms but excludes the $uu_{xx}$ and $u_{x}^{2}$ terms at $\sigma=0.01\text{ and }0.05$. For larger values of $\sigma$, no correct terms are included in the final equation form.

We investigated if the recovered terms from the PDE-FIND with pruning algorithm using ANN approximations can be used as the specified mathematical model in an inverse problem methodology (cf., [7]) to recover the final parameter estimate values from the nonlinear Fisher-KPP Equation. If we take the union of all terms that are included in the final equations in Supplementary Table S3 for the ANN method using noisy data ($\sigma>0$), then we have an equation of the form

[TABLE]

We estimated the parameters $a,...,f$ in Equation (18) for each value of $\sigma$ by simulating the solution to this PDE using the method of lines and minimizing Equation (8) using the Nelder-Mead algorithm. We input the equations from Supplementary Table S3 for the ANN method as the initial guess for each data set. We find that performing this inverse problem leads to accurate parameter estimates for the true terms in the nonlinear Fisher-KPP Equation and small coefficient values for the incorrect terms ($u_{xx}$ and $1$) for $\sigma=0,0.01,0.05,\text{ and }0.25$ (Table 4). At $\sigma=0.10$ and $0.50$, this inverse problem methodology leads to small coefficient estimates for $uu_{xx}$ in addition to $u_{xx}$ and $1$. Note that this same process would not lead to ultimately recovering the true equation and parameter estimates from the spline or finite difference approximations because their final equations never included the correct $uu_{xx}$ term in the final equation for noisy data ($\sigma>0$).

4 Conclusions and Future Work

The novel use of the ANN method presented here is a significant step toward making PDE learning more achievable in realistic scenarios with noisy biological data. Because an ANN is a fully differentiable function, it can be used to approximate derivative computations to build the library of terms needed for PDE learning. The current practice to build a library of terms for learning PDEs when noise is present in the observed data $u(x,t)$ is to use finite difference or spline approximations for small amounts of noise [28, 35]. Our findings suggest that these two methods are highly sensitive to the amount of noise in the data in the range of less than $5\%$ noise. Moreover, these methods can not incorporate heteroscedastic statistical models such as Equation (2) to handle non-constant error noise, which is a typical phenomenon encountered in biological data. We showed that the ANN method outperforms spline and finite difference approximations of $u(x,t)$ and its partial derivatives when significant levels of non-constant error noise are present in the data $u(x,t)$.

In general, we found that when using the ANN method, the true underlying equation is more accurately recovered than when using finite differences or splines at all noise levels considered in this work (up to $\sigma=0.50$). The disparity between the accuracy of the ANN and the other two methods sharply increased with the noise levels in the data. For the diffusion-advection equation, the ANN allows PDE-FIND to learn the correct equation for up to 25% noise levels, while splines and finite difference calculations fail for noise levels over 10%. The ANN data allows PDE-FIND to learn the Fisher-KPP Equation for up to 10% noise levels whereas the other two methods fail at 5%. All methods struggle with the nonlinear Fisher-KPP Equation, but the ANN data leads to an equation that can then be used with an inverse problem methodology [7] to infer which terms are meaningful.

There are several reasonable explanations for the difficulties in learning the nonlinear Fisher-KPP Equation. Three of the terms in Equation (6) are the product of two terms including $u$ or its derivatives ($u^{2}$,$u_{x}^{2}$, and $uu_{xx}$). In practice, these terms may be inaccurate approximations from noisy data (as demonstrated in Table 3). Multiplying two inaccurate terms may lead to an even larger amount of uncertainty associated with these estimates. We postulate that the high level of uncertainty in these type of terms resulting from the product of inaccurate estimates likely increases the difficulty of learning to include them in the process of PDE learning. Furthermore, it must be noted that the data for this equation was generated numerically the by finite difference method. Though some analytical solutions to the Fisher-KPP equations are known, they come either in series form, which introduces error by series truncation, or in traveling wave form, which would be indistinguishable from the advection equation to PDE-Find. The finite difference methods used to approximate the spatial derivatives in the nonlinear Fisher-KPP equation introduce second and fourth-order truncation error terms which lead to numerical dispersion and diffusion effects that may account for the recovery of some unexpected terms.

We found that the use of pruning following our implementation of the PDE-FIND algorithm increased our ability to recover the correct equation in terms of the TPR. It may be argued that these additional terms learned from the PDE-FIND implementation without pruning would be removed from the final equation if more regularization (i.e., a larger value of $k$ in implementation of the Greedy algorithm) were used. However, when faced with the issue of learning the governing equation from actual data in practice, one will not have the ability to know when the specified equation is correct or not. We thus need to identify the correct hyperparameters without any a priori knowledge. We observed that the systematic biases in our ANN (depicted in Figure 2) make it difficult to choose a hyperparameter value that leads to the correct equation because these biases are present in the training and validation data. The pruning algorithm is a way to correct for when the incorrect hyperparameter has been chosen by ensuring that all terms in the learned PDE are sufficiently sensitive to constitute a strong signal in the data, instead of resulting from a bias in our approximation methods.

Our use of pruning to remove terms from learned PDEs could be improved in future work. While effective, our implementation is somewhat crude, in which we pre-specify a threshold level to prune parameters based on out-of-sample MSE values on the validation data set. Previous studies have discussed F-statistics as one way to infer the increase in variance that pruning a variable will lead to, but there are still many different interpretations of these results which makes a definitive statistical pruning method challenging to ascertain [11].

In future work, we may also consider extending the use of the ANN method to statistical models that a relevant to other types of biological noise, such as weighted least squares statistical models [7] or log-normally distributed errors [9]. We will also investigate using our methodology in the case where the statistical model is not known and needs to be estimated from data with a model-free numerical method [4].