Global optimization using Gaussian Processes to estimate biological parameters from image data
Diana Barac Michael D. Multerer Dagmar Iber

An Example on s-H-Convexity in ℂ^2
    Berit Stensønes
    December 30, 2023
==================================



Keywords: parameter estimation; global optimization; Gaussian processes; image-based modeling; parameter estimation for computationally costly models

Declarations of interest: None.





§ ABSTRACT

Parameter estimation is a major challenge in computational modeling of biological processes. This is especially the case in image-based modeling where the inherently quantitative output of the model is measured against image data, which is typically noisy and non-quantitative. In addition, these models can have a high computational cost, limiting the number of feasible simulations, and therefore rendering most traditional parameter estimation methods unsuitable. In this paper, we present a pipeline that uses Gaussian process learning to estimate biological parameters from noisy, non-quantitative image data when the model has a high computational cost. This approach is first successfully tested on a parametric function with the goal of retrieving the original parameters. We then apply it to estimating parameters in a biological setting by fitting artificial in-situ hybridization (ISH) data of the developing murine limb bud. We expect that this method will be of use in a variety of modeling scenarios where quantitative data is missing and the use of standard parameter estimation approaches in biological modeling is prohibited by the computational cost of the model.





§ INTRODUCTION

Computational models of biological systems typically serve two purposes: to test the validity of underlying assumptions and to make predictions about the behaviour of the system <cit.>.  Such models often contain parameters whose true values are unknown. This can be because the parameters correspond to physical constants e.g. production rates, that are difficult or expensive to measure accurately, or it can be that the parameter has no direct biological equivalent that can be determined experimentally. The behaviour of the model, and hence the predictions derived from it, can differ depending on the parameter set used. Consequently, parameter estimation is an important problem in computational modeling, see <cit.> and the references therein. 

Recent advances in imaging technologies, which enable ever increasing spatial and temporal resolution, have resulted in widespread use of image data in computational modeling <cit.>.  The unique advantage of image data such as in situ hybridization (ISH) or immunohistochemistry (IHC) data lies in is its ability to reveal the spatial and temporal distribution of a molecule of interest in an anatomically intact system <cit.>. However, in most cases, this data is non-quantitative and noisy, making its integration into a quantitative computational model challenging. In addition, biological processes are generally comprised of a large number of coupled components that react in a non-linear fashion and whose behavior may vary in both space and time. This spatio-temporal behaviour can be modelled by coupled non-linear partial differential equations (PDEs), e.g. reaction-diffusion equations, whose solution usually entails a high computational cost. In this paper, we address the problem of parameter estimation for image-based modeling in which the model is computationally expensive and the data is noisy and non-quantitative. 

Parameter estimation aims to find the set of parameters for which the model best fits the available data. To this end, the first step consists of defining a suitable cost function which measures the distance between the model output and the available data, and is therefore dependent on the type of data used. When quantitative data is available, a common choice for the cost function is the sum of squared deviations (least squares) <cit.>. We look at the case where the model output and available data are images, thus, it is necessary to define an appropriate image metric that compares the non-quantitative data with the quantitative model output. The image Euclidean distance (applying the sum of squared deviations to image pixels) is the most commonly used image metric <cit.>. We describe a way to define the cost function, taking into account the non-quantitative nature of the image data. To this end we weaken this metric to a pseudometric by decomposing the image into areas of high and low intensity, prioritizing that the model recapitulates those rather than a set gradient, which can vary between images. 

The next step in parameter estimation is the minimization of the cost function with respect to the parameters. If the cost function has a single minimum (or the parameters can be restricted to a range such that only one minimum is present) local optimization methods can be used. These can be gradient-based or direct, see e.g. <cit.> and the references therein. However, as computational models are often complex and nonlinear, this generally leads even simple cost functions, such as least squares, to possess multiple minima. Then, it is necessary to use a global parameter estimation method as a local one will likely result in the algorithm getting trapped in a local minimum. The necessity of global estimation methods in biological modelling is discussed in <cit.>.

There exist several algorithms for global optimization, see <cit.> and the references therein. The earliest and simplest of these methods is multi-start (running a local method repeatedly from a set of predefined initial starting points). However, this approach usually does not work for realistic applications, due to its computational inefficiency. Simulated Annealing is widely applied in global optimization problems and was shown in  <cit.> to give the best solution in a comparison of parameter estimation methods to estimate rate constants. Nonetheless, the authors noted the high computational cost of the algorithm. Evolutionary Computation methods (including Genetic Algorithms, Evolutionary Programming and Evolution Strategies) are widely used in global optimization <cit.>. Based on the concept of genetic evolution, these methods take an initial population of parameters and iteratively refine it (usually by selecting the best parameters and recombining them in some way) until an optimum parameter set is found. Hence, these methods also require a large number of individual simulations. Hybrid methods, i.e. the combination of different optimization methods, have also been used successfully in parameter estimation of biological models <cit.>.

An alternative approach is to infer the posterior distribution over the parameters (the probability of the parameters given the data) instead of just retrieving a point estimate. This can be achieved using Bayesian methods <cit.>. As the calculation of the likelihood function is often analytically and computationally intractable, approximate Bayesian computation (ABC) methods, which replace this calculation with a comparison between observed and simulated data, are a popular choice. ABC methods based on Markov chain Monte Carlo (MCMC) and sequential Monte Carlo (SMC) have been successfully applied to biological systems <cit.>. As these involve sampling a large number of parameters to estimate the posterior distribution, these approaches are usually only feasible for lower-dimensional problems where model simulation is computationally cheap.

As can be seen, most of the global optimization methods used in biological parameter estimation require a large number of independent model simulation steps, making them well suited for models that are computationally cheap to solve, but impractical for models whose simulation is computationally expensive. Surrogate or response surface modeling is an approach to global parameter estimation that is applied when the number of possible model evaluations is greatly limited <cit.>. Response surface methods aim to estimate the response surface (the surface of the cost function as a function of the parameters) in an adaptive fashion and usually consist of two steps that are iterated <cit.>. In the first step, the response surface is estimated using the model outputs that have been computed thus far. The second step consists of selecting a suitable next point for evaluation. This is usually based on a tradeoff between exploration (selecting a point in a region where the model has not been evaluated) and expected improvement (selecting a point that is most likely to give a lower value of the response surface). Different functions can be used to estimate the shape of the response surface from the given data. One such class is regression models, where the response surface is assumed to take a particular shape (e.g. quadratic or linear) and the surface is fitted accordingly. This works well if the assumed shape is indeed an accurate assumption of the response surface, but poorly otherwise <cit.>. Another class is interpolating models, where the response surface is made up of basis functions based on the value of the response surface at known parameter values. Gaussian processes are one example of such a method <cit.> and can be used to estimate the response surface and the uncertainty of that estimate at each point, thereby facilitating step 1 and step 2 of the response surface method. Although Gaussian processes are applied in biological modeling, their use is rather restricted to predicting the state variables of a model <cit.>, than directly estimating the cost function with respect to the parameters. In this paper, we apply Gaussian processes to estimate the cost function and thereby estimate the parameters of our model. 

Parameter dependencies are frequently present in biological models <cit.>. If a model parameter cannot be unambiguously estimated from the available data it is said to be non-identifiable. Determining the presence of parameter dependencies is important as model predictions may differ depending on the model parameters i.e., two sets of parameters that fit the available data equally well may give rise to different model predictions. Non-identifiable parameters can be either structurally or practically non-identifiable. A parameter is structurally non-identifiable if the model structure prevents it from being estimated unambiguously, independent of the quality of the data. A parameter is practically non-identifiable if the amount or quality of the available data (rather than the structure of the model) prevents its unambiguous estimation. Structural non-identifiability is reviewed in <cit.> and both structural and practical non-identifiability are discussed in <cit.>.

For biological systems described by ordinary differential equation (ODE) models, various methods for detecting structural non-identifiability have been proposed <cit.>. For arbitrary models, the profile likelihood approach, a general approach to identify structural and practical non-identifiability, is introduced in <cit.>. The identifiability of the parameter is investigated by varying it and looking at the minimum objective function value at each point under all possible values of the remaining parameters. However, if the model is not analytically tractable and is computationally costly, calculating this profile likelihood may not be feasible. We use multivariate sensitivity analyses to determine parameter dependencies, and therefore parameter identifiability, in our models as the above approaches are not appropriate due to the models' “black box" nature and high computational cost.

The rest of the paper is structured as follows: in the next section we describe the parameter estimation process. This is split into two parts: the first considers defining the cost function in the absence of non-quantitative data and the second describes the parameter estimation process. The following section presents the results. To test the method, we first use it with images generated from a parametric function with the goal of retrieving the original parameters. We then present the result of estimating parameters in a more complex biological setting with the aim of fitting artificial in-situ hybridization (ISH) data of the developing murine limb bud. The final section is the discussion where we consider the impact of non-identifiability, the potential limitations of the process and future work. 



§ METHOD



 §.§ Cost Function

The first step of the parameter estimation process consists of defining the cost function. We develop a method for the case where the model is solved on a domain similar to to that of the shape of the tissue of interest from which our data is obtained, where the data to be reproduced is non-quantitative, noisy image-data. An example of such data is ISH data. ISH is a technique in which a complementary strand to the nucleic acid of interest, with a reporter molecule attached, is applied to the tissue of interest <cit.>. Upon addition of a substrate, a chemical reaction catalysed by the reporter molecule takes place which reveals the location of a specific segment of nucleic acid, for instance DNA or RNA. Variations in the experimental protocol, for example leaving the probe to bind for a different amount of time, washing for a longer period, or, in particular, the length of time of the chemical reaction, result in differences in the data obtained. In addition, non-specific binding of the probe leads to background signal. The result of this is non-quantitative and noisy image data, wherein the relationship between the intensity of the staining and the concentration of the nucleic acid of interest is non-linear, see Figure <ref>A. 

From the image data we determine a reference image, the shape of which is identical to the domain the model is solved on. We identify three separate areas in the data: the area where the staining is consistently darkest in all the data sets, i.e. the concentration is highest, the area where staining is consistently low across all data sets and the remaining area, where the staining is neither high nor low or varies between data sets. These three areas are marked on the reference image (the high concentration area in black, the low in white and the rest in gray), see Figure <ref>B. As the data is noisy and non-quantitative the gradient of gene expression can vary between images of the same gene (Figure <ref>A). However, usually it is clear which areas of the tissue of interest have reproducibly higher levels of gene expression and which have lower levels of expression. Therefore, we want our model to recapitulate those areas of high and low gene expression, rather than a set gradient, which can vary between images. 


We then quantify the difference between the output of the model and the reference image. The reference image is grayscale i.e. black and white. White corresponds to a maximum pixel value of 255/255, and black corresponds to a minimum pixel value of 0/255. Shades of gray take on pixel values within this range. In order to compare the output of our model to the reference image we normalize the output of each component in the model we wish to compare with respect to its maximum concentration and plot the result in grayscale, with black (pixel value 0/255) representing the maximum concentration and white (pixel value 255/255) zero concentration. We now have two grayscale images that can be compared. 

We define the cost function as follows. Areas corresponding to the white area in the reference image are penalized linearly if their pixel value is smaller than 179/255 i.e. if the concentration is more than 30% of the maximal value. This means that the penalty is zero if the pixel value, v, is within the range of 179 to 255 and otherwise equal to 179-v. Similarly areas corresponding to the black area in the reference image are penalized linearly if their pixel value is larger than 76/255 i.e. the concentration is less than 70% of the maximal value. Areas corresponding to the gray areas in the reference image are not included in the calculation. The cost function is then the squared sum of errors across all pixels across all images (one image for each variable we wish to include in the cost function), see Listing <ref>.

By splitting the reference image into areas of high, mid and low concentration and normalizing with respect to the maximum concentration we remove the impact of noise from the data on the cost function and the problem of comparing non-quantitative data to the quantitative output of the model. Other thresholds are, of course, possible as well.

[frame=single, caption = Cost Function, label=1st:label, language = C, basicstyle = , mathescape=true]  

for images i
	for pixel indices p in image i
		if referencePicture(p,i)≤76
			/* 70
			scorePic(p,i) = max(testPicture(p,i)-76,0)
		else if referencePicture(i,j)≥179
			/* 30
			scorePic(p,i) = max(255-testPicture(p,i)-76,0)
		else
			scorePic(p,i) = 0
		end
	end
end
/* score = sum of squared entries of scorePic */
score = ∑_i,pscorePic(p,i)^2








 §.§ Parameter Estimation

We use a Gaussian process to estimate the response surface with respect to our parameters. A Gaussian process is defined as a collection of possibly infinite random variables such that any finite collection of them has a multivariate normal distribution. A Gaussian process is fully defined by its mean function, μ(x,θ) and covariance kernel function k(x,x';θ).  In this case x, x'∈ X, where X is the parameter space. The kernel parameters θ are usually dependent on a characteristic length scale. The characteristic length scale defines roughly how far away two points have to be in order to be almost uncorrelated, see <cit.> for an overview of Gaussian processes. With the mean and covariance function defined, one can then calculate the expected value of the response surface, which we call f, at a set point, given the value of f at other points. It is also possible to calculate the standard deviation of f at that point, which tells us how likely it is that the value of f at that point will be close to the expected value calculated. 

We use the  algorithm in MATLAB <cit.> for our parameter estimation process which makes use of Gaussian processes to predict the response surface. It uses the automatic relevance determination (ARD) Matern 5/2 kernel. This is defined as

    k(x,x' | θ) = σ_f^2 (1+ √(5) r + 5/3r^2 ) exp(-√(5) r),

where

    r = √(∑_m=1^d(x_m-x_m')^2/σ_m^2).

Here, θ is the vector of the kernel parameters. In this case

    θ_m    = logσ_m, for m = 1, 2, ...d  and
    θ_d+1   = logσ_f,

where d is the number of parameters that are being estimated, σ_f is the standard deviation of the response surface and theσ_m are the characteristic length scales for each parameter.
 estimates the parameters θ of the covariance function from known values of the response surface. It then uses Gaussian process regression to predict the shape of the response surface and selects which point to sample at next. The theoretical convergence properties of this approach are analysed in <cit.> and <cit.>.

As the model to be solved is assumed to be computationally expensive, there is a limit to the number of feasible iterations. We therefore set a limit to the number of iterations we allow for the parameter estimation procedure. This number will differ depending on the model e.g. a model with a high number of parameters to be estimated or a highly irregular response surface will generally require more iterations for the identification of a suitable parameter set. A model that has a very high computational cost will have a lower number of feasible iterations than a model whose computational cost is less prohibitive. The set number of iterations is a trade-off between the accuracy of the parameter estimation process and the computational cost of the process. 

For our model we limit the number of iterations to 400. Estimating a response surface is based on the knowledge of the values of the cost function at the parameter values sampled thus far. Therefore, we first run the model at parameter values spaced evenly in the sample space to gain a rough estimate of the shape of the response surface. We use a Halton sequence <cit.> to generate 100 points that are spaced evenly in the parameter space. To generate these points we use the MATLAB function . We choose to use a Halton sequence rather than selecting points randomly as a Halton sequence distributes points consistently in space (unlike a sequence of points generated randomly) <cit.>, see Figure <ref>. See <cit.> for details on generating Halton sequences. We input the parameter values at which we have already evaluated the model with the corresponding value of the cost function into the  function as initial (already known) values. We allow the algorithm to run for a further 300 iterations (resulting in 400 iterations in total) and then select the parameter with the lowest corresponding cost function value that has been sampled by the algorithm as our optimal parameter.





§ RESULTS




 §.§ Estimating the parameters of a test function



To test our parameter estimation procedure we used it to  estimate the coefficients of the function

    f( x,θ)=max{0,1-c| x_2-p( x,θ)|},   c>0,

where

    p(α,θ) = θ_0+∑_k=1^M-1
    +θ_2kcos(kα)+θ_2k-1sin(kα)+θ_2M-1cos(Mα),     0 ≤α < 2π

is a trigonometric polynomial and α = arg x. We set the parameter c to 4. Thus, the parameter estimation amounts in this case to estimating the 2M coefficients θ of p(α,θ)
from images. We set M=4 and fix θ_2=0 (excluding it from the parameter estimation procedure). We then try to estimate the remaining 7 parameters ρ = [θ_0, θ_1, θ_3, θ_4, θ_5, θ_6, θ_7]. We ran the parameter optimization procedure with a total of 400 iterations, of which the first 100 were initial Halon points. The range of the parameter space we sampled across was [0.01,1] for each parameter with the Halton set distributed across the log of the parameter space. This means that, for all parameters, the number of Halton points present in the range [0.01,0.1] should be approximately equal to the number present in the range [0.1,1]. The parameter estimation procedure was also performed on the log of the parameter space as we assumed changing a parameter from 0.01 to 0.1 would have a similar effect to that of changing a parameter from 0.1 to 1. We found that for 7 parameters we were able to successfully reconstruct the image, with the algorithm returning an optimum parameter set close to the original one, see Table <ref>. We then tested the algorithm on a more blurred image by changing the value of c to 2 (corresponding to more ambiguous or noisy image data) and the algorithm was still able to reconstruct the image. See Figure <ref>. To examine the effect of parameters that have little impact on the cost function we added dummy parameters into the algorithm that had no effect on the images. These were input into  as parameters to optimize but were not inputs to the parametric function. The algorithm was still able to find a value close to the original parameter set for the 7 relevant parameters even with 10 dummy parameters added, see bottom row of Table <ref>.








We then looked at how decreasing the number of total iterations and changing the number of initial Halton points (seed number) affected the optimization procedure. We ran the optimization procedure with 200, 300 and 400 total iterations, with 0, 25, 50 and 100 initial Halton points. We found that simulations run with no initial Halton points performed considerable worse and those with 200 and 300 iterations did not perform as well as those run with 400 iterations. The results are shown in Figure <ref>A.

We further examined how varying the number of Halton points and the number of points  used to predict the response surface affected the optimization procedure for a total of 400 iterations. The active set is defined as the set of points used to predict the response surface. The size of this set can be modified by changing the variable GPActiveSetSize in . For example, if  has sampled at 350 points but the active set size is 300 then 300 points out of the 350 are selected uniformly at random without replacement to estimate the response surface. (The previous simulations had been run with the default GPActiveSetSize=300.) We tested each of these points with an active set size of 300 and an active set size of 400. The results are shown in Figure <ref>B. We found that the lowest median value was for 100 Halton points with GPActiveSetSize = 400. However, simulations run with 25 and 50 Halton points also performed well, whereas those with 200 and 0 Halton points performed considerably worse. This indicates that, for a set number of iterations, an initial set of Halton points can greatly improve the parameter estimation procedure, but too many prevent the algorithm from finding the optimum.

Finally, we looked at how varying the number of Halton points and the number of iterations performed after these initial simulations affected the parameter estimation procedure. We call the latter the number of bayesopt iterations. The total number of iterations is then the sum of the number of initial Halton points and the number of bayesopt iterations. We varied the seed number (number of initial Halton points) in increments of 100, from 0 to 300. We varied the number of bayesopt iterations from 100 to 400, also in increments of 100. The results are shown in Figure <ref>C. As expected, with the number of Halton points fixed, the parameter estimation procedure generally performed better the higher the number of bayesopt iterations. Unexpectedly, the converse does not appear to be true. With a fixed number of bayesopt iterations, increasing the seed number does not lead to a smaller minimum objective value.  For the values sampled, the optimum seed number for a fixed number of iterations was 100. This optimum seed number likely depends on the number of parameters to be estimated, the range of the parameters and the number of bayesopt iterations and should ideally be determined before the parameter estimation procedure. This could be done by running the procedure on a surrogate model which has the same number of parameters to be estimated as the original problem but is cheaper to compute.




 §.§ Comparison with other algorithms

We then tested how this procedure compared to two other commonly used optimization algorithms: simulated annealing and  genetic algorithm. The functions  and  in MATLAB use simulated annealing and genetic algorithm respectively to search for a function minimum. We ran both algorithms on our test function to try and retrieve the original parameters with 400, 700 and 1000 iterations. The results are shown in Table <ref>. As can be seen, neither  nor  succeeded in retrieving the original parameters of the test function, even when the number of iterations was increased to 1000. 


max width=,keepaspectratio

algorithm     iterations     time ()    θ_0     θ_1     θ_3     θ_4     θ_5     θ_6     θ_7     objective 
 N/A     N/A     N/A     0.5     0.2     0.1     0.3     0.3    0.2     0.1
 
     400     2215 ± 515    0.5 ± 0.01     0.19 ±0.02     0.08 ± 0.02     0.30 ± 0.02     0.30 ± 0.02     0.19 ± 0.01    0.08 ± 0.02     14293003

    400     602 ± 17     0.32 ±  0.12     0.55 ± 0.23     0.38 ± 0.23     0.57 ± 0.19     0.38 ± 0.19     0.46 ± 0.2     0.28 ± 0.17     516660821

    700     1047 ± 90     0.35 ± 0.13     0.54 ± 0.22     0.31 ± 0.19     0.56 ± 0.18     0.37 ± 0.2     0.44 ± 0.21     0.24 ± 0.17     390425288

    1000     1448 ± 5     0.35 ±  0.13     0.50 ± 0.23     0.17 ± 0.1     0.53 ± 0.18     0.23 ± 0.17     0.28 ± 0.12    0.11 ± 0.12     254855944

    400     725 ± 13     0.55 ± 0.05     0.1 ± 0.11     0.02 ± 0.02     0.26 ± 0.08     0.32 ± 0.04     0.09 ± 0.02     0.08 ±  0.05     564173347

    700     1152 ± 11     0.55 ±  0.03     0.11 ± 0.02     0.05 ± 0.03    0.25 ± 0.04     0.31 ± 0.04     0.13 ± 0.07    0.07 ± 0.03     405229775

    1000     1615 ± 31     0.55 ± 0.06     0.13 ± 0.07     0.05 ± 0.05     0.26 ± 0.04     0.32 ± 0.03     0.12 ± 0.05     0.08 ± 0.04     383413571




Comparison of global parameter estimation algorithms. Original parameters of the parametric function (top row) and optimum parameters found by different algorithms (sample size = 10). As the algorithms are not deterministic, distinct optimal values can emerge in independent runs. Objective is average minimum objective function value found. Simulation time of parametric function = 1.42 ± 0.05s (sample size = 10).




 §.§ Estimating Network Parameters




We then tested our procedure on a more complex, biologically relevant model with 7 parameters. We based our model on the developing murine limb and fitted it to artificial ISH data of Fgf10 expression. During early limb development Fgf8 is expressed at the apical ectodermal ridge (AER), a structure at the distal end of the limb bud. It diffuses into the mesenchyme where it binds to its receptor to form a signaling complex (F8R). This complex upregulates the expression of another protein (Fgf10), see Figure <ref> for a network scheme and the corresponding ISH expression data. We solved our model on a realistic growing domain with domain size (initial size:160  from AER to flank and final size: 240  from AER to flank) and time (5 hours) corresponding to the development of the limb bud from somite stage 21 to somite stage 24. Fgf10 expression at the time point corresponding to somite stage 24 was used to compare against Fgf10 expression evaluated at the same time point from artificial ISH data (generated with parameter set θ_0).

The following reaction-advection-diffusion equations are used to represent the gene expression dynamics with the terms to be interpreted as follows. The variables C_F8, C_F10 and C_F8R represent the concentration of FGF8, FGF10 and the FGF8-receptor-ligand signalling complex respectively. As the receptor is only present in the mesenchyme <cit.>, we define the signalling complex C_F8R to be present exclusively in this part of the computational domain. FGF10 and FGF8 are diffusible and can therefore be present in both the mesenchyme and AER domain. p_F10 and p_F8 are the production rates of FGF10 and FGF8, respectively. We set our initial values to zero. The initial expression of C_F10 is then also zero. As Fgf10 is expressed early in the limb bud, we add an additional production term to C_F10: p_0Fgf10. The indicator functions 1_Mes and 1_AER are used to show that Fgf10 and Fgf8 are expressed in the mesenchyme and AER respectively. The diffusion rate of FGF10, D, is assumed to be equal to that of FGF8 as they are both diffusible ligands. The diffusion rate of the signalling complex is denoted D_r. Similarly, d is the degradation rate of FGF8 and FGF10, and d_r that of the signalling complex. FGF8 binds to its receptor at a rate of k_on8 and the unbinding rate of the receptor-ligand complex is k_off8. As the FGF8 receptor is only present in the mesenchyme, binding and unbinding only occur in that part of the domain. The total receptor concentration is assumed to be in quasi-steady-state, we therefore use RT_F8 to denote the receptor concentration, so (RT_F8-C_F8R) is the concentration of free receptor. Assuming mass action kinetics, the binding and unbinding terms are 1_Mesk_on8C_F8(RT_F8-C_F8R) and 1_Mesk_off8C_F8R respectively. Finally, we model C_F8R-induced expression by a Hill function C_F8R^2/C_F8R^2+K_F8R_ F10^2 with Hill constant K_F8R_ F10^2 and Hill coefficient 2. The velocity field of the tissue, due to growth, is denoted u. 


    ∂ C_F10/∂ t + ∇(uC_F10)    = D  C_F10  + 1_Mes(p_F10C_F8R^2/C_F8R^2+K_F8R_ F10^2+ p_0Fgf10)-d C_F10,
    
    ∂ C_F8/∂ t + ∇(uC_F8)    = D  C_F8  + 1_AERp_F8- d C_F8-  1_Mesk_on8C_F8(RT_F8-C_F8R) +1_Mesk_off8C_F8R,
    
    ∂ C_F8R/∂ t + ∇(uC_F8R)    =  D_r C_F8R  + k_on8C_F8(RT_F8-C_F8R)-(k_off8+d_r)C_F8R.


While the length and time scales of the developmental processes are typically well known, the absolute value of protein concentrations are largely unknown (due to ISH being non-quantitative). The standard procedure would be to non-dimensionalize the model to reduce the number of parameters. However, given the available information regarding the developmental length and time scales, we will retain dimensional length and time scales, and only eliminate two parameters by scaling the variables with respect to their particular production rates, i.e.


    C_F10=p_F10C_F10,  C_F8=p_F8C_F8, C_F8R=p_F8C_F8R,  and RT_F8=p_F8RT_F8R.



We get 


    ∂C_F10/∂ t + ∇(uC_F10)     = D C_F10  + 1_Mes(C_F8R^2/C_F8R^2+K_F8R_ F10^2/p_F8^2+  p_0Fgf10/p_F10)- dC_F10,
    
    ∂C_F8/∂ t + ∇(uC_F8)    = D C_F8  +1_AER - dC_F8 +  1_Mes(k_off8C_F8R -k_on8p_F8C_F8(RT_F8-C_F8R)), 
    
    ∂C_F8R/∂ t + ∇(uC_F8R)    =  D C_F8R  +k_on8p_F8C_F8(RT_F8-C_F8R) -(k_off8+d_r)C_F8R.


We assume k_off8≪ d_r,  k_on8p_F8. Then, grouping parameters by settting: K_F8R_ F10^2 = K_F8R_ F10^2/p_F8^2,  k_on8= k_on8p_F8 and p_0Fgf10=p_0Fgf10/p_F10 and dropping the tildes for convenience we obtain the final set of equations:


    ∂C_F10/∂ t + ∇(uC_F10)     = D C_F10  + 1_Mes(C_F8R^2/C_F8R^2+ K_F8R_ F10^2+p_0Fgf10)- dC_F10,
    
    ∂C_F8/∂ t + ∇(uC_F8)     = D C_F8  +1_AER- dC_F8+1_Mes( -k_on8C_F8(RT_F8-C_F8R)), 
    
    ∂C_F8R/∂ t + ∇(uC_F8R)    =  D_r C_F8R  +k_on8C_F8(RT_F8-C_F8R) -d_rC_F8R.


The diffusion rate of the ligand, D, is fixed as 100 times higher than the diffusion rate of the receptor D_r. The artificial ISH data of Fgf10 corresponds to the production rate term of FGF10, i.e. 1_Mes(F8R^2/F8R^2+ K_F8R_ F10^2). The 7 network parameters to be estimated are then as follows: the diffusion coefficient for FGF8 (D), the ligand degradation rate (d), the signaling complex degradation rate (d_r), the binding rate between FGF8 and its receptor (k_on8), the total amount of receptor (RT_F8) and the Hill coefficient for enhanced production of Fgf10 (K_F8R_ F10^2). The seventh parameter is the initial FGF10 production rate (p_0Fgf10). We took an initial parameter set θ_0 (see Table <ref>) and used it to generate a reference image for expression of Fgf10 at the time point corresponding to somite stage 24: see Figure <ref>. We then ran our parameter estimation process, using the in silico data to calculate the cost function. We ran the parameter estimation procedure with 100 Halton points and 400 iterations in total with an active set size of 400. The optimum parameter set found by the algorithm resulted in an expression pattern that was almost identical to the input image, see Figure <ref>. However, the optimum parameter value that was returned varied significantly from our initial parameter set θ_0, with only the degradation rate of the signaling complex within 10% of the original parameter value, see Table <ref>.


Comparison of original parameters (top) of the network and optimum parameters found (bottom)

D=θ_1 (^2^-1)    d=θ_2 (^-1)     d_r=θ_3 (^-1)    k_on8=θ_4 (^-2)    RT_F8=θ_5 ()    K_F8R_ F10=θ_6 ()    p_0FGF10=θ_7 
 16000     0.1     0.3     60     5     0.15     0.1
 
 16344     0.216     0.627     87.7     2.38     0.101    0.0146





In order to understand why the algorithm failed to retrieve the original parameter values, although the image representing Fgf10 expression at the optimum value found was almost identical to the image representing Fgf10 expression for the initial parameter set θ_0, we ran sensitivity analyses for the model. We also ran sensitivity analyses for the parametric function, where the parameter estimation process had accurately retrieved the original image as well as the original parameter set.  

The results of the single variable sensitivity analysis can be seen in Figure <ref>. The response curve for the parametric function is convex about the optimum and continuous. For the parameters in the limb network the response curve for θ_2 is discontinuous. For θ_7 it is flat, implying that this parameter has no impact on the cost function. For parameters θ_5, θ_4 and θ_1 the response curve appears almost flat for a wide subset of parameters about or on one side of the optimum. This could be a reason as to why the optimum parameters found for parameters θ_4, θ_7, θ_2 and θ_5 are different to the initial θ_0. 

The results of the multivariate sensitivity analysis, for selected parameters, can be seen in Figure <ref>B. The impact of the parameters of the parametric function on the cost function is independent. On the other hand, parameter dependencies can be seen for the limb network model. For example, if we look at parameters θ_3 and θ_5 (d_r  and RT_F8) we see that if θ_3 is fixed at a higher value, the minimum value of the cost function is obtained by decreasing θ_5. If we look at the optimum parameter set that was returned and the original values of θ_3 and θ_5 we can see this. For the optimum value returned, θ_3 is higher than the original value and θ_5 is correspondingly lower. Similarly, for parameters θ_3 and θ_6 (d_r  and  K_F8R_F10), for the optimum value returned, θ_3 is higher than the original value and θ_6 is correspondingly lower.

We further looked at whether it was possible to simplify the model. To this end, we implemented a direct activation of C_F10 by C_F8, thereby removing C_F8R from the model and simplifying the system. This resulted in the following equations:


    ∂ C_F10/∂ t  + ∇(uC_F10)    = D  C_F10  + 1_Mes(C_F8^2/C_F8^2+ K^2 + p_0Fgf10)- dC_F10,
    
    ∂ C_F8/∂ t  + ∇(uC_F8)    = D  C_F8  +1_AER- dC_F8


with the 4 parameters: D,  d,  K and  p_0FGF10. We then ran the parameter estimation procedure, using the same reference image as for the original model, and performed a multivariate sensitivity analysis about the optimum found. The optimum parameters found are shown in table <ref> and the sensitivity analysis in Figure <ref>C. Interestingly, the optimum parameters found were close to the corresponding original parameters of the original network. However, the image of Fgf10 expression was considerably different to that of the reference image, see Figure <ref>D. This is likely due to the difference in diffusion rate between F8 and F8R, leading to a sharper boundary in the reference image, resulting in a qualitatively different result when F8R is removed from the equations. From the sensitivity analysis it can be seen that parameter dependencies are still present, even in this simplified model, for instance between D and d and D and K.



  








§ DISCUSSION

The importance of global optimization methods for parameter estimation in biological modelling is well recognized (due to the complexity of the cost function with respect to the parameters) and numerous methods have been developed for this purpose. Most require a large number of model simulations to obtain good results. This is not feasible if the model is computationally costly, thereby prohibiting multiple iterations. In this paper, we proposed a heuristic global parameter estimation procedure for image-based modelling where the image data is noisy and non-quantitative and the model is computationally expensive to solve. To this end, we described a cost function, which measures the distance between the current and desired output of the model, and a procedure, using Gaussian process learning, to subsequently perform the parameter estimation. 

Using an image-generating parametric function to test the process we showed that, for seven parameters, this method retrieves the original parameters (and therefore the original image) when the parameters are independent. We then tested the procedure on a biological network with artificial data. As is often the case in biological models <cit.>, even with simple networks, there was a high level of parameter dependence in the network model. Applying the parameter estimation procedure resulted in a parameter set corresponding to an output image closely resembling the reference image. However, due to the parameter dependence, the optimal parameter set found differed considerably from that of the parameter set used to generate the reference image. 

There are several factors which influenced the parameter dependence in our model. For dynamical models, such as the limb model in our paper, a necessary condition to fix parameter information is the availability of temporal information. If the time scale is unknown, there will necessarily be parameter dependence as parameter sets that result in the same behaviour of the model, differing only in the rate of that behaviour, will be equally valid. In this case, parameter dependence due to lack of temporal information can be removed by non-dimensionalizing the equations with respect to time. In our limb model, we only evaluate the similarity of the model to the data at one time point (somite stage 24). This is one reason we have high parameter dependence: as parameter sets that fit the data equally well at this point, but result in different gene expression patterns at later points, have the same cost function value. Evaluating at additional time points, for instance at later somite stages, would decrease the parameter dependencies in our model, albeit at an increased computational cost due to increased simulation time of the model. In particular, this should decrease the dependencies between the degradation of the receptor d_r and the diffusion coefficient D and the degradation of the receptor d_r and the Hill constant K. 

In addition, the only read-out we have is the production rate of FGF10 which we compare against artificial expression data. Including read-outs for other variables would decrease the parameter dependence, as well. For example, there is parameter dependence between K and parameters that control the concentration of F8R e.g. d_r,  k_on8,  RT_F8 as the production rate of FGF10 can be modified by either changing the sensitivity of FGF10 to F8R (via K) or by changing the concentration of F8R. Both the sensitivity K and the concentration of the FGF8-receptor complex are difficult to measure. But if data was available such that, for instance, the concentration of F8R could be added to the cost function, then the parameter dependence between K and RT_F8 and K and d_r should decrease. The level of parameter dependence is contingent on the amount of data that is available to evaluate the goodness of fit of the model. In cases where the concentration of all the components in the model is known at several time points, the parameter dependence should be low. However, this is rarely the case in biological modelling, where either experiments can be expensive to run (prohibiting measurements at several time points) or substances can be difficult to measure (in which case the concentration of several components in the model will be unknown). 

Other than testing the validity of underlying assumptions of the biological system, the main purpose of computational modeling is to obtain predictions. Although additional experiments can be performed to constrain model parameters <cit.>, in most cases, the degree of measurement accuracy required to constrain predictions is extremely high and therefore impractical <cit.>. Therefore, the impact of parameter dependence on model behaviour should be examined. If the model behaviour is insensitive to certain parameter combinations, then parameter dependence does not rule out accurate predictions. In this case, the model output that determines model predictions is as insensitive to certain parameter combinations as the model output that is used to fit the model, and predictions can be made despite parameter non-identifiability <cit.>. Model predictions under different sets of feasible parameters can be analysed (the ensemble method) <cit.>. For example, uncertainty estimates for the model output can be calculated by generating an ensemble of feasible parameters and calculating the mean and standard deviation of model output. This approach is illustrated in <cit.>, where experimental verification confirms model predictions, despite large parameter uncertainties in the model. However, in cases where the model is computationally expensive, determining and simulating the sets of feasible parameters could have a prohibitive computational cost. 

When considering which global optimization method to use, an important consideration is the computational cost of the model. The simulation time of our limb model was approximately 40 seconds, running 400 simulations therefore took approximately 4.5 hours. The Gaussian process regression algorithm has a time complexity of O(n^3), where n is the number of training data. Hence, if the simulation is relatively cheap, the ratio of computational time of the algorithm to the computational time of the simulations can be high. In this case, global optimization algorithms with lower computational complexity, such as genetic algorithms, can perform better. However, as mentioned in the introduction, these algorithms often require a high number of iterations to find an optimum. 

We used the function  in the MATLAB optimization toolbox to test how genetic algorithms would perform on our parametric function. When the function was run with a maximum of 400 iterations the optimum found by  was far from the original input values. The image obtained with the optimum found was also significantly different from the original image. With 700 iterations  again failed to retrieve the original parameters, however, the image obtained by the optimum found was similar to the original image. Using the function  with 700 iterations of our parametric function took, on average, 1152 seconds. As one model simulation, on average, of our parametric function was 1.42s this meant that 994s could be accounted for from just running the simulations. Therefore,  spent approximately 158s on operations other than running the simulation (a relatively small fraction of the total running time). In contrast, running  with 400 iterations of our parametric function took 2215s. This meant that 568s could be accounted for from just running the simulations, and then  spent approximately 1647s on operations other than running the simulation (a considerably larger fraction of the total running time). Consequently, for our parametric function, running 700 iterations of  is less costly than 400 iterations of . In contrast, the simulation time for the limb model is 40 seconds. Running  with 400 iterations then takes approximately 400 × 40s + 1647s=17,647s. In comparison, running  with 700 iterations would take a minimum of 700 × 40s + 158s=28,158s. Accordingly, the computational cost of the model determines which global optimization methods are feasible to employ.

In this paper, we restricted our models to those with seven parameters or less. To apply this method to estimate a larger number of parameters most likely a larger number of iterations would be required, and would need to be determined, along with an appropriate ratio of Halton points (initial seeds) to total iterations. The total number of iterations required could become prohibitively large for models with a large number of parameters and high computational complexity. A solution to this could be having a large number of initial seeds which are computed in parallel (if this is possible), thereby reducing the number of iterations that have to be performed in serial for the rest of the optimization procedure. However, as can be seen in Figure <ref>C, having a larger number of initial seeds does not necessarily lead to better results. It would therefore be useful to determine the optimum seed number (for the number of parameters to be estimated and the number of feasible bayesopt iterations) and then run these initial seeds in parallel. 

To calculate our cost function, we generated a reference image from the image data available. This reference image discriminates areas of high and low concentration, as well as the area that should be excluded from the cost function calculation. Depending on the features of the data that the model output should match, a different number of areas can be identified and penalised accordingly. For example, if a subset of the image should have a medium intensity, this area can also be identified and incorporated into the cost function, by penalising model outputs that result in concentrations that are too high or too low. What percentage of the maximum concentration is considered a high or low concentration depends on the model and data and the cost function can be implemented accordingly. The parameter estimation procedure can also be applied in cases where a portion, or all, of the data is quantitative. In this case, the cost function is simply modified by penalising the difference between the model output and the absolute concentrations expected rather than the maximum concentration output of the model. 

We expect that the method described in this paper can be used to perform global parameter estimation when model simulation is computationally costly and quantitative data is missing. The step of generating the reference image from the available data could be automised using image processing. This would facilitate an unbiased determination of the cost function from the data and the reference image could be updated as new data becomes available. The method could also be used as a model selection tool <cit.>, by treating model structure as a parameter, where different network topologies are compared to see which produces the best fit to the data.



§ SOFTWARE

The parametric function was solved using MATLAB (version R2017b)  <cit.>. The model of the developing limb was solved with COMSOL Multiphysics 5.1 <cit.> and MATLAB Livelink (version R2017b) <cit.> as described in <cit.>.



§ ACKNOWLEDGMENTS

We thank Harold Gómez and Marcelo Boareto for critically reading the manuscript, and Guido Sanguinetti for discussions. This work was supported by a SystemsX.ch iPhD grant (2014/240 to D.I.) and the Swiss National Science Foundation (SNF.170930 to D.I.).


unsrt

