Integration of in vitro and in silico Models Using Bayesian Optimization With an Application to Stochastic Modeling of Mesenchymal 3D Cell Migration

Cellular migration plays a crucial role in many aspects of life and development. In this paper, we propose a computational model of 3D migration that is solved by means of the tau-leaping algorithm and whose parameters have been calibrated using Bayesian optimization. Our main focus is two-fold: to optimize the numerical performance of the mechano-chemical model as well as to automate the calibration process of in silico models using Bayesian optimization. The presented mechano-chemical model allows us to simulate the stochastic behavior of our chemically reacting system in combination with mechanical constraints due to the surrounding collagen-based matrix. This numerical model has been used to simulate fibroblast migration. Moreover, we have performed in vitro analysis of migrating fibroblasts embedded in 3D collagen-based fibrous matrices (2 mg/ml). These in vitro experiments have been performed with the main objective of calibrating our model. Nine model parameters have been calibrated testing 300 different parametrizations using a completely automatic approach. Two competing evaluation metrics based on the Bhattacharyya coefficient have been defined in order to fit the model parameters. These metrics evaluate how accurately the in silico model is replicating in vitro measurements regarding the two main variables quantified in the experimental data (number of protrusions and the length of the longest protrusion). The selection of an optimal parametrization is based on the balance between the defined evaluation metrics. Results show how the calibrated model is able to predict the main features observed in the in vitro experiments.

Introduction

Cell migration is a fundamental event in a wide variety of physiological processes, spanning from embryogenesis (Knecht and Bronner-Fraser, 2002; Martin and Parkhurst, 2004), angiogenesis (Lamalice et al., 2007; Spill et al., 2015), osteogenesis (Reina-Romo et al., 2012), inflammatory response (Luster et al., 2005), immune response (Bogle and Dunbar, 2010), and wound healing (Shaw and Martin, 2009; Valero et al., 2014a), to develop diseases such as cancer and metastasis formation (Franz et al., 2002; Condeelis et al., 2005; Condeelis and Pollard, 2006).

Cell migration can present different characteristics according to the dimensionality in which it is produced. Thus, cell migration on 2D surfaces has been widely studied and is typically characterized by a balance between counteracting traction and adhesion forces (Sunyer et al., 2016; Escribano et al., 2018). However, cells generally migrate in a 3D extracellular matrix (ECM) adopting different migration strategies regulated by several factors such as the cell type and the properties of the ECM. In these 3D environments, the mechanisms governing cell migration are far less understood due to both the technical challenges and the complexity of migratory behaviors (Zhu and Mogilner, 2016).

Based on the cell type and the cellular microenvironment (Te Boekhorst et al., 2016; Talkenberger et al., 2017)—in particular ECM parameters such as density, porosity and stiffness—, individual cells migrate using two distinct mechanisms (Friedl and Wolf, 2010; Swaney et al., 2010; Bear and Haugh, 2014). When cells are unable to adhere to the ECM, they modify their shape and squeeze through the ECM pores by using the amoeboid migration, which is very efficient—rapid cell locomotion (cell speed ~10 μm/min)—and it is observed in cells such as neutrophils and T cells (immune system) (Beauchemin et al., 2007; Lämmermann et al., 2008; Swaney et al., 2010). In contrast, whenever cells adhesion to the ECM is high, they degrade the ECM to pass through by using the mesenchymal migration mode, which is very inefficient—cell displacement is very slow (cell speed < 1 μm/min)—and it is observed in cells such as fibroblast (wound healing) and osteoblasts (bone formation) (Friedl and Wolf, 2010). This mesenchymal migration mode is investigated in this paper.

In vitro experiments have become increasingly sophisticated in order to reproduce as accurate as possible the natural biological surroundings of organisms from in vivo studies. As in vitro studies have increased their sophistication, their requirements have also grown in complexity. Due to the complexity and the expensive lab work of in vitro experiments, in silico studies have a complementary role in understanding mesenchymal cell migration. Computer-based mathematical models allow performing a vast number of controlled and reproducible experiments with much lower associated costs. In fact, these computational models can be classified according to several factors such as the numerical approach of the biological processes: continuous (Vermolen and Javierre, 2012; Valero et al., 2014b; Serrano-Alcalde et al., 2017), discrete (Meineke et al., 2001; Bentley et al., 2009; Scianna et al., 2012; Scianna and Preziosi, 2014; Van Liedekerke et al., 2015), or hybrid (Alber et al., 2007; Bauer et al., 2009; Daub and Merks, 2013; Milde et al., 2014; González-Valverde and García-Aznar, 2018). In addition, they can also be classified according to the type of movement that cells develop as individual (Schlüter et al., 2012; Trichet et al., 2012; Ribeiro et al., 2017; Moure and Gomez, 2018), if cells migrate independently, or collective, forming an interconnected meshwork or cluster (Bazmara et al., 2015; González-Valverde et al., 2016; Norton and Popel, 2016; Camley and Rappel, 2017; Escribano et al., 2018). Computational models can also be classified as mechanical (Zaman et al., 2005; Borau et al., 2011), biochemical (Hatakeyama et al., 2003; Provenzano et al., 2008), or mechano-chemical (Kim et al., 2015, 2018; Moure and Gomez, 2017; Ribeiro et al., 2017).

More recently, different authors (Sun and Zaman, 2017; Kim et al., 2018; Mark et al., 2018) have focused their works in the combination of in vitro experiments and in silico modeling in order to elucidate the influence of specific factors on individual and collective cell migration. The combination of both methodologies opens new opportunities for research, because models allow the simulation of in vitro conditions in order to directly obtain additional information not available from experiments, but that can be indirectly evaluated in-vitro. For example, recently, Sunyer et al. (2016) analyzed collective cell durotaxis, combining experiments with numerical models in order to understand that the difference of stiffness sensed by cells at both edges of the cell monolayer promotes the directional migration.

In this work, we propose to establish a new strategy based on the Bayesian optimization (BO) technique, which combines numerical simulations relied on a mathematical model and in vitro experiments in order to calibrate the model's parameters. In particular, a mechano-chemical model of individual mesenchymal 3D migration is presented, with a focus on accelerating the numerical simulations that determine the 3D migration trajectories. This strategy allows the full integration of numerical models and experimental measurements in order to improve knowledge of how cells regulate this mesenchymal 3D migration.

Materials and Methods

This section is organized in order to describe how experimental measurements and numerical simulations can be integrated in a consistent way. To facilitate their explanation, first, we briefly describe the mathematical model of cell migration (Ribeiro et al., 2017) and its numerical implementation. Next, we show the results from in-vitro experiments and their quantification. Then, we present how both results can be combined by means of Bayesian optimization in order to calibrate the numerical model with the experimental results. Finally, we test the potential of our calibrated numerical model under different chemoattractant concentrations and gradients.

Model Description

The selected model to simulate 3D cell migration is based on a previous one (Ribeiro et al., 2017) (Figure 1). Here, we describe the main aspects of this model in order to understand how the full calibration of this model is developed. This model assumes that cell migration can be described by three clearly differentiated stages. During the first stage the cellular chemosensing mechanism allows cells to probe the chemical cues located on their surroundings through different membrane receptors (Roca-Cusachs et al., 2013; Moreno-Arotzena et al., 2015). In particular, the focus is on how fibroblasts detect molecules of the chemo-attractant factor Platelet-derived Growth Factor (PDGF from now on) through a specific cell surface receptor, the tyrosine kinases one (also known as RTK) (Cao et al., 2004; Poukkula et al., 2011). The second stage simulates how the activation of these receptors triggers intracellular processes that regulate the onset of dendritic protrusions in different directions throughout the ECM (Weiger et al., 2010; Liou et al., 2014). In fact, these protrusions can protrude (pushing the matrix) and contract (pulling the matrix). Lastly, the third stage models how the dynamics of these protrusions regulate cell migration in 3D (Campellone and Welch, 2010; Starke et al., 2013; Moreno-Arotzena et al., 2014) by establishing a relation between the contractile force generated by each protrusion and the cell body translocation.

FIGURE 1

Figure 1. Global model scheme. (Left) the chemosensing mechanism simulates how RTKs located in the cell membrane become activated by binding to PDGF-BB molecules (blue circles). This RTK activation, in turn, triggers the activation of PI3K molecules inside the cell (PI3K inactivated molecules as gray triangles and PI3K activated ones as red triangles). (Middle) protrusions (${\text{p}}_n^i$) grow and stabilize on those areas with high concentration of PI3K activated molecules. (Right) the longest protrusion generates a traction force (${\text{F}}_{trac}^{lp}$) when retracting, which exerts a reaction force (${\text{F}}_R^{lp}$) over the cell body. As a result of these reaction forces, the ECM generates a drag force (F_drag) over the cell body.

Next, these three main stages of the process of cell migration are described in greater detail. But first, the model of 3D cell behavior is defined.

Modeling Cell Behavior

The 3D structure of the cell is geometrically modeled as a set of one-dimensional bars representing dendritic protrusions (Ribeiro et al., 2017). Those bars are located in a three-dimensional environment and diverge from a central connecting point that represents the cell body. This central connecting point—which can be associated to the cell nucleus or the cell centrosome—exists solely for modeling purposes as the point where all the bars are connected (Figure 1 right).

Modeling the Chemosensing Mechanism

This first stage models how the chemically reacting system that allows the cell to sense the chemo-attractant factor (located in the surrounding ECM) evolves through time (Figure 1 left).

It is assumed that the only signal pathway guiding protrusion dynamics is the one including just a chemo-attractant factor located in the ECM, RTKs in the cell membrane and PI3K molecules inside the cellular body. The PDGF has been chosen as the chemical factor to interact with the cell due to its pivotal role in regenerative processes (Chen et al., 2007; Friedlaender et al., 2013; Elangovan et al., 2014; Shah et al., 2014). However, the model could be extrapolated to other growth factors.

In order to replicate the cellular chemosensing mechanism, our model simulates the interaction between different species through time. From a temporal perspective, the simplified mathematical model that mimics this chemosensing mechanism is based on a set of reactions (Equation 1) (Hatakeyama et al., 2003; Hawkins et al., 2006) and defined by a set of differential equations (Equation S1 of the Supplementary Material).

$\begin{array}{cc}\{\begin{array}{c}RTK+{\left[F\right]}_{k_2}{\rightleftharpoons }^{k_1}RTKF\text{ }(R_1 R_2)\\ RTKF+PI3K{\to }_{k_3}RTKF+PI3K_A\left(R_3\right)\\ PI3K_A{\to }_{k_4}PI3K(R_4)\end{array}& (1)\end{array}$

From a spatial perspective, it is assumed that membrane receptors such as RTKs are homogeneously distributed over the cell surface. However, the activation density of such membrane receptors depends on the distribution of chemoattractant molecules (F). In particular, there are more activated RTK receptors (RTKF) on those areas of the cell surface surrounded by a higher concentration [F]. In contrast, on those areas of the cell surface surrounded by a lower concentration [F], there are less activated receptors. Thus, cells are able to sense the spatial distribution of F.

Two sources of stochasticity in cell migration are associated to the chemosensing mechanism: the evolution of the chemically reacting system (defined by Equation 1 for the proposed model) and the activation of RTKs based on the concentration of chemoattractant molecules surrounding the cell. Therefore, the chemical reactions defined by Equation (1) are assumed to be stochastic processes described by a Poisson distribution (Ueda and Shibata, 2007). This premise makes possible to consider receptor activation over a domain with varying concentration of factor F to be a multivariate non-homogeneous Poisson's distribution. Therefore, it is possible to model this activation of RTKs by means of the Inverse Method described by Saltzman et al. (2012).

By computing an approximate solution of this problem, it is possible to evaluate at any given time (t_k) the variation of PI3K_A in any specific location of the cell surface. In order to estimate this spatio-temporal variation we defined the variable s(α, β, t_k) which stores the spatial persistence of PI3K_A activation across time (t_k), in a space location of the cell surface defined by coordinates (α, β) (Figure 2)—since we are dealing with a 3D model of a cell, we represent the cell membrane as a flat surface defined by the polar coordinates α and β. Therefore, the signal s = s(α, β, t_k) is evaluated at a fixed time t_k by means of the convolution function, taking into account the temporal evolution of the chemical signal in this surface location and its surroundings—roughly an area the size of a protrusion section. The model equations guiding the chemosensing mechanism are summarized in Table 1.

FIGURE 2

Figure 2. (Left) spherical coordinates (α, β) of point q. (Right) s signal distribution example over the cell membrane.

TABLE 1

Table 1. Equations associated with the chemosensing mechanism.

We assume that cell's consumption of chemoattractants is negligible. Thus, the chemoattractant chemical profile does not change with time.

Modeling Protrusion Dynamics

Once the tempo-spatial variation of activated PI3K (PI3K_A, s) is estimated, it is possible to determine protrusions location by means of a set of thresholds (s_birth, s_exp, and s_ret) that act as a signal filter. In particular, s_birth represents the minimal amount of signal s that cells need to develop new protrusions, as suggested by many authors (Ueda and Shibata, 2007; Weiger et al., 2010; Jilkine and Edelstein-Keshet, 2011; Chen et al., 2017); those points inside the cellular body where s is higher than s_birth are considered locations where novel protrusions sprout. Furthermore, any pre-existing protrusion becomes reinforced if, in its location, s is higher than s_exp. However, if there is not enough signal s for the protrusion to remain active, it becomes unstable; in those points where s is lower than s_ret pre-existing protrusions retract and disappear (Table 2). Besides, in order to simplify the search of signal s peaks where protrusions centroids are localized, an internal model parameter (s_binary) is used to transform s into a binary signal. This means that only during the protrusions localization, any surface point where s is lower than s_binary becomes 0 whereas every surface point with s greater or equal to s_binary becomes 1.

TABLE 2

Table 2. Equations associated with the protrusions dynamics.

In addition, it is assumed that this signal variation δs also regulates, in conjunction with the ECM mechanical properties, the protrusive stretch characteristics due to the cytoskeleton activity.

Therefore, protrusions generate forces against the ECM. Consequently, the mechanical properties of the ECM act as a regulator for the extension or retraction of protrusions, as suggested by Liou et al. (2014) (Figure 3). This behavior is simulated by considering protrusions analogous to an elastic inclusion (ellipsoid) embedded in the ECM, applying Eshelby's analytical solution of ellipsoidal elastic inclusions in an elastic infinite body (Eshelby, 1957). We consider that during this second stage protrusions grow inside a collagen-based fibrous matrix and they adhere to ECM fibers. Thus, we consider the ECM behaves as a linear elastic material that constrains the growth of protrusions. In fact, during this growth, protrusions push to the ECM deforming it and the elastic properties of the ECM regulate this deformation. In this case we quantify the growth of the protrusion and the deformation of the ECM by means of the Eshelby's theory, assuming the protrusion as an inclusion that is embedded in the ECM. Moreover, in all the cases we assume infinitesimal deformation.

FIGURE 3

Figure 3. Scheme of the three different configurations in protrusion dynamics. ${\epsilon }_k^f$ represents the free expansion/retraction (ECM does not restrict protrusions deformation) Cauchy's strain tensor. ε_o^* is the compatibility Cauchy's strain tensor. ${\epsilon }_k^c$ represents the total deformation Cauchy's strain tensor. We assume infinitesimal deformation.

Equations guiding protrusions dynamics are summarized in Table 2.

Modeling Cell Body Translocation

Finally, based on the experimental observations of how protrusions determine cell body translocation (Moreno-Arotzena et al., 2015; Del Amo et al., 2017; Movilla et al., 2018), it is assumed that the longest protrusion determines cell motion directly. The longest protrusion presents a larger adhesion surface and, consequently, adhesion proteins have higher probability to connect with the ECM. Every cell protrusion, except the leading one, becomes non-adherent and, as a result, they are all dragged by the cell during cell motion. The retraction of the leading protrusion generates a reaction force ($F_R^{lp}$) supported by the cell body. Thus, by focusing just on the reaction force generated by the longest protrusion ($F_R^{lp}$), it is possible to estimate the exerted drag force (F_drag) by the ECM on the cell body (Figure 1 right). As a result, both cell speed and position can be estimated at any given time t following the definition proposed by Borau et al. (2011), which takes into account the ECM viscosity. During the third stage we model the cell body translocation and, as the position of the cell center is modified, we assume the cell body is on the fluid component of the ECM. Thus, we consider that the cell is moving through a fluid. As a result, and in order to compute the drag force exerted by the ECM on the cell body, we take into account the viscosity of the ECM. The equations guiding cell body translocation are summarized in Table 3.

TABLE 3

Table 3. Equations associated with the body translocation.

We assume that there is a mechanical balance between the traction force of the adherent protrusion ($F_{trac}^{lp}$), the longest one, and its corresponding reaction force ($F_R^{lp}$) supported by the cell body due to $F_{trac}^{lp}$ (Figure 1 right). Equation (12) defines a relationship between the contractile force magnitude ($F_{trac}^{lp}$), due to actomyosin activity, and the protrusion length.

Numerical Implementation

Our computational model has been designed using a scheme based on the three fundamental mechanisms: chemosensing mechanism, protrusions dynamics, and the cell body translocation (Figure 4). These three stages have been implemented in Python using powerful packages and libraries for scientific computing such as Numpy (van der Walt et al., 2011) and SciPy (Jones et al., 2001) to maximize the model's performance.

FIGURE 4

Figure 4. Global algorithm scheme. First, the chemosensing mechanism is simulated using the tau leaping algorithm. In this first stage, the concentration and gradient of the PDGF-BB is the main influence factor. During the second stage of the process, it is taken into account both ECM mechanical properties and cell mechanics in order to simulate protrusions development. Finally, the modeling of the cell body translocation is also influenced by the ECM mechanical properties (in particular, ECM viscosity) as well as cell mechanics. The blue boxes are associated with the chemical factor, the red ones with cell mechanics and the yellow ones with the ECM mechanical properties.

The stochastic time evolution of the given set of reactions (R₁, R₂, R₃, and R₄) had been numerically simulated by using, originally, the Stochastic Simulation Algorithm (SSA; also known as the Gillespie Algorithm) (Gillespie, 1976, 1977) in the first version of this work (Ribeiro et al., 2017). However, the SSA is considered too slow for our purposes and a faster alternative is proposed, the tau-leaping algorithm (Gillespie, 2001; Cao et al., 2006). The SSA computes an exact solution of the time evolution of a chemically reacting system. In contrast, the tau-leaping algorithm estimates a good enough¹ approximation (Lok, 2004; Cazzaniga et al., 2006) by leaping over many reactions at once using Poisson random numbers.

The tau-leaping method tries to accelerate stochastic simulations by approximating the frequency of each reaction being fired in the next specified time interval [t, t+τ]. By comparison, the SSA focuses only on one reaction per time interval which may be prohibitively small (Anderson et al., 2011). As long as the value of τ is small enough so the leap condition² is satisfied, it is possible to compute a good approximation of the evolution of a given chemically reacting system.

It is worth to mention that neither the SSA nor the tau-leaping algorithm use a fixed time step to simulate the evolution of biologically reacting systems like the one presented in this work. Instead, they compute a new value τ in each iteration based on the current state of the system and a random variable.

The initial amounts of each reactant as well as the reaction rates (k₁, k₂, k₃, and k₄) used are included in Table 4.

TABLE 4

Table 4. Initial amounts of each reactant as well as the reaction rates obtained from literature.

Based on the spatial distribution of PI3K_A molecules as well as their concentration on those locations, protrusion growth is then set. Protrusion final length is computed by applying Eshelby's solution of ellipsoidal elastic inclusions in an infinite elastic body. Mechanical equations are analytically solved using a computational algorithm. An elastic modulus of 104 Pa is assumed for the ECM based on previous experimental works of gels with a concentration of 2 mg/ml collagen type I (Movilla et al., 2018; Valero et al., 2018).

Lastly, the mechanical equilibrium associated to protrusion-generated forces is solved. Then, taking into account that the longest protrusion is the one leading cell migration, it is computed both cell speed and position in the following time increment.

We decouple the simulation of the chemosensing mechanism from the other two stages of the model (protrusions dynamics and cell body translocation) because we are considering two different time scales in our model. In fact, the chemical and mechanical events occur at different time scales. In order to accurately simulate the proposed chemically reacting system we are using the iterative tau leaping algorithm with a variable associated time step τ in the range [0.5, 1.5] seconds. However, to model protrusion dynamics and the cell body translocation we are using a time step dt of 5 min. Actually, it is because variations in the signal s (Equation 2) between two consecutive time steps t and t + τ are really small whereas protrusions require more noticeable variations of the chemical signal in order to change their current state. As a result, it is required to keep track of the cumulative variation δs.

Development and Quantification of in vitro Experiments

Once we have numerically implemented the proposed model, we have to calibrate its parameters in order to optimize the performance of this computational model. We calibrate the model here presented by comparing the results of its simulations with experimental data. In particular, we focus on two different features to fit the model's parameters: the length of the longest protrusion and the number of protrusions of the migrating cell. As a result, we have performed in vitro studies to get accurate experimental measurements of the length of the longest protrusion and the number of protrusions.

In vitro experiments have been performed by culturing Normal Human Dermal Fibroblasts (NHDF)—human skin primary cells—within 2 mg/ml collagen gels at a concentration of 2.5x10⁵ cells/ml, and temperature and atmosphere conditions have been maintained at 37°C and 5% CO₂. Immediately after the seeding, cells' evolution has been monitored with multidimensional microscopy for 4 h (from 0 to 4 h), every 5 min and 5 μm of Z axis, with 200x magnification (20x objective) and phase contrast (Figure 5). We have chosen a 2 mg/ml collagen concentration because it already implies a matrix pore size (< 1 μm) (Fraley et al., 2010). Individual cell protrusions have been quantified by in-house Matlab algorithms (Moreno-Arotzena et al., 2015). For each image stack, best Z has been chosen in order to maximize accuracy and minimize the complexity of the manual analysis of both the cell body and its protrusions. Single cell analysis of four different samples has been performed for the given collagen concentration (2 mg/ml).

FIGURE 5

Figure 5. Norman Human Dermal Fibroblast (NHDF) cultured in 3D collagen-based fibrous matrix (2 mg/ml collagen). Image was captured with a Nikon D-Eclipse Microscope with a Plan Fluor 200x magnification (20x Objective) and phase contrast.

FGM™-2 (Fibroblast Growth Medium-2) has been used to support the growth of primary human fibroblasts. It contains a supplementation of GA-1000, recombinant human insulin 0.5%, HFGF-B GF, and 2% of Fetal Bovine Serum. Thus, these in vitro experiments only include a very low and fixed concentration of growth factors included in the culture medium; they do not include any chemoattractant gradient.

Model Calibration Using Bayesian Optimization

During the last couple of decades, as the available computational power has greatly increased, so has the complexity of in silico models and the number of parameters included in those models. As a result, the complexity of the calibration process has also increased. However, it is still often the case that this calibration process is performed using a very manual approach. Each parameter must be tuned manually despite the search space being usually too vast to be effectively navigated. Besides, there may be interactions or dependencies between some parameters. This process can be very tedious, especially when dealing with in silico simulations that require several hours of execution time.

This calibration process can be mapped to a non-linear optimization problem where the objective is to find the simulation parameters that best fits the in vitro experiments. In this way, we are able to automate the process. However, most non-linear optimization solvers require a large number of iterations, gradient information of the fitting function or they are sensible to local optima. In our case, the large number of iterations could make the problem intractable as the evaluation of the fitting function associated to our in silico model is very costly because it requires several simulations of our stochastic model.

More formally, we are looking for the set of optimal experiment parameters x that satisfy:

$\begin{array}{cc}\begin{array}{r}{x}^{\star }={\displaystyle \underset{x\in \chi }{\text{argmax}}}f(x),\end{array}& (14)\end{array}$

where f is the fitting function between the in vitro and the in silico models and X is the parameter search space as defined in Table 5.

TABLE 5

Table 5. Model parameters calibrated using Bayesian optimization with SigOpt.

Bayesian optimization, also called Efficient Global Optimization (EGO) (Jones et al., 1998) is a general purpose black-box optimization methodology that it is characterized for requiring a very small number of iterations before reaching global optimization. Thus, it is especially suitable for experimental design and calibration of expensive processes (Shahriari et al., 2016). Bayesian optimization uses a probabilistic surrogate model of the target function combined with optimal decision theory to drive the search toward the global optimum in less iterations than popular non-linear optimization alternatives like PSO (Kennedy and Eberhart, 1995), CMA-ES (Hansen et al., 2003) or L-BFGS (Nocedal, 1980). In the case of Bayesian optimization, the surrogate model uses machine learning to capture previous iterations acting as a memory of the full optimization process. Meanwhile, the decision component carefully selects the next query at each iteration.

In the case of simulation calibration, there are many variables that can be used for fitting, some of them might be even competing. Then, we can redefine the problem as one of multi-objective, multicriteria optimization or Pareto optimization:

$\begin{array}{cc}{x}^{\star }=\underset{{}_{x\in \chi }}{\arg \max }(f_1(x),f_2(x),\dots ,f_n(x)),& (15)\end{array}$

In this case, the objective is not only to find a single optimal value for the simulator parameters, but to find the whole set of Pareto optimal points, that is, those points that dominate the rest of the points. Although this is a completely different problem, the seminal work of Knowles (2006) extended the Bayesian optimization methodology to the multi-objective setup.

There are several pieces of software that implements Bayesian optimization, like BayesOpt (Martinez-Cantin, 2014). A full review can be found in Shahriari et al. (2016). However, many of them do not support multi-objective optimization and those that do support multiple criteria are very limited in terms of other features. In this work, we have used SigOpt³ (Martinez-Cantin et al., 2018) for its support for parallelization and multi-objective optimization. Besides, it provides other features like the parameter importance, which will be discussed in the Results section.

For our experiments, we have decided to fit two competing metrics: the length of the longest protrusion (llp) as well as the number of protrusions (np) (Figure 1). The fitting of the in silico values with respect to the in vitro measurements is computed using the Bhattacharyya coefficient (also known as BC), which has been widely used to compare the similarity or discriminate of two continuous or discrete distributions (Comaniciu et al., 2000). In fact, for discrimination it corresponds to the upper bound of the Bayesian error when performing Bayesian hypothesis testing with symmetric cost functions and uninformative priors (Nielsen, 2014). Note that, Bayesian hypothesis testing already includes a penalization for model complexity and priors result in a regularization effect, being less sensitive to overfitting than classical hypothesis testing (Kass and Raftery, 1995).

In particular, histograms of both in vitro and in silico experiments are used as discrete distributions to compute those metrics (Equation 16).

$\begin{array}{cc}\begin{array}{r}BC={\displaystyle \sum _{i=1}^n}\sqrt{his{t}_{invitro}^i·his{t}_{insilico}^i},\end{array}& (16)\end{array}$

histⁱ represents the value of the i-th histogram bin defined as the probability of occurrences of values in the range (x_i−1, x_i].

The selection of metrics affects model calibration, so we have carefully selected the metrics with a greater influence on cell migration to the best of our knowledge. Moreover, these metrics are based on experimental measurements that we are able to accurately quantify. However, there are other measurements based on cell motion, such as the instant cell speeds, that are so low that we are not able to quantify them with the required accuracy. For those metrics it is only possible to perform a qualitative analysis. Although our proposed metrics are based on just two quantities measured in the experimental data: the length of the longest protrusion and the number of protrusions, we consider that both variables are fundamental in the regulation of the final 3D cell motion. In particular, experimental observations (Moreno-Arotzena et al., 2015; Del Amo et al., 2017; Movilla et al., 2018) suggest that the length of the longest protrusion has great influence over the cell speed whereas the number of protrusions has a great impact on the cell trajectory (whether it is random or directional).

Optimizing the BC function can be considered as a form of Bayesian learning in the sense that we are trying to fit a model that best represents the distribution of the data, and therefore maximizing the posterior. Similarly, optimizing the BC can be seen as a form of Bayesian hypothesis testing where we are rejecting all the models with higher Bayesian error.

Furthermore, Bayesian optimization is a black-box method, meaning that it does not require specific knowledge about the metric, and that metrics can be easily interchanged. Thus, the same methodology can be applied to any other feature or any other similarity metric, such as KL-divergence or any other loss function. Besides, we can include metrics not directly related to the data such as cost, time, etc. These metrics can be competing, meaning that one metric cannot be improved without another metric suffering. As a result, the solutions distributed in the Pareto set might be distributed in a complex way. Thus, sample efficient search like Bayesian optimization is of paramount importance. Besides, the resulting Pareto front allows the expert user to balance the competing metrics a posteriori, choosing the most convenient parametrization in different circumstances.

Model Validation Using Different Chemoattractant Concentrations and Gradients

After calibrating the numerical model, we have to validate it, testing their predictive ability to simulate different cell responses under different chemical gradients. This validation process allows us to prove that the proposed model does not only accurately replicate the results used to calibrate it, but also new ones, so that there has been no overfitting during the calibration process. In the preceding calibration process, we have used quantitative results related to both the length of the longest protrusion and the number of protrusions of migrating NHDF from in vitro experiments without any chemoattractant gradient. However, the validation process of this computational model is based on qualitative observations of migrating cells surrounded by a chemoattractant factor diffusing throughout the ECM (Song et al., 2006; Bosgraaf and Van Haastert, 2009). We have simulated six different extracellular environments. Three of these extracellular environments include different PDGF gradients (10⁻¹, 10⁰, 10¹ μM/mm) but a fixed PDGF concentration at the initial cell's position (0.8 μM). The other three extracellular environments include a fixed PDGF gradient (10⁰ μM/mm) but different PDGF concentrations at the initial cell's position (0.08, 0.8, and 8.0 μM). Twenty simulations have been executed for each extracellular environment, using the same seeds used during the calibrating process. The comparison between in vitro and in silico results is based on qualitative observations of the velocity component in the direction of the chemotactic gradient (v_x).

We assume a fixed growth factor profile without any induced modifications of the spatial gradient due to the growth factor diffusion throughout the ECM. Thus, the chemoattractant chemical profile is assumed to be temporally stable as the inlets and outlets of our system keep a fixed growth factor profile during our 4-h simulation.

Results

By means of in vitro experiments in fibroblasts, it is possible to quantify both the length of every protrusion, as well as the number of protrusions generated at every checkpoint t (t = 0, 5, 10, …, 240 min). Figure 6 shows an example of the images generated by multidimensional microscopy and the posterior protrusions analysis performed using in-house Matlab algorithms. However, our model focuses on the length of only the longest protrusion at each temporal checkpoint t, ignoring the length of the other protrusions, as explained in Section Modeling Cell Body Translocation. Therefore, during the calibration process the comparison between in vitro and in silico experiments is performed by means of the BC using these two features (length of the longest protrusion and the number of protrusions generated by migrating cells).

FIGURE 6

Figure 6. (Left) Phase contrast example of a NHDF cell cultivated in a 2 mg/ml collagen gel, with 200x magnification (20x objective) acquired using multidimensional microscopy. (Right) Protrusion analysis performed by in-house Matlab algorithms; red line delimits cell body, yellow lines represent the protrusions, and blue line shows cell body displacement. In this case, the longest protrusion is the green one and the number of protrusions is 5.

During calibration, for every iteration in the optimization loop, 20 simulations replicating the in vitro scenario of a 2 mg/ml collagen ECM have been executed—in order to capture the stochastic nature of our model. Those 20 simulations used 20 different seeds in order to initialize the global random number generator of our model. Once the 20 simulations have been completed, their associated histograms are computed by means of a computer-based algorithm. These histograms (e.g., Figure 7 bottom) are compared with the in vitro histograms (Figure 7 top) using the proposed evaluation metrics BC_llp and BC_np (defined in Equations 18, 19 respectively and based on Equation 16).

$\begin{array}{cc}\begin{array}{ccc}B{C}_{llp}& =& \frac{{{\displaystyle \sum }}_{i=1}^{N}B{C}_{llp}^i}{N},N=20,\end{array}& (18)\end{array}$

$\begin{array}{cc}\begin{array}{ccc}B{C}_{np}& =& \frac{{{\displaystyle \sum }}_{i=1}^{N}B{C}_{np}^i}{N},N=20,\end{array}& (19)\end{array}$

In order to compute the two metrics using the BC, it is required to generate the associated histograms for both the longest protrusion length and the total number of protrusions. Histograms associated to in vitro experiments using a cellular microenvironment based on 2 mg/ml collagen gels show how the protrusion length ranges from over 0 μm to almost 140 μm. However, most of the longest protrusions have a length in the interval 40–60 μm (Figure 7 top left). Regarding the number of protrusions, there is a high dispersion, ranging from 1 to 14 protrusions in each individual fibroblast during migration (Figure 7 top right). Figure 7 (bottom) shows an example of a couple of histograms associated to the in silico experiments. In this case, we have generated in silico histograms using the best parametrization suggested by SigOpt with metrics BC_llp = 0.87 and BC_np = 0.81 (Figure 7 bottom). These histograms show how, although the length of the longest protrusions is between 0 and more than 150 μm, there is a peak in the interval 60–80 μm (Figure 7 bottom left). Regarding the number of protrusions, there are usually about 9 to 12 in each fibroblast during migration (Figure 7 bottom right). When comparing measurements of the length of the longest protrusion, the mean values are 63.71 (in vitro) vs. 65.98 (in silico), whereas the standard deviations are 31.20 (in vitro) vs. 26.82 (in silico). For the measurements of the number of protrusions, the mean values are 7.57 (in vitro) vs. 7.38 (in silico), whereas the standard deviations are 3.27 (in vitro) vs. 4.00 (in silico).

FIGURE 7

Figure 7. Normalized histograms associated to in vitro experiments (Top) based on the length of the longest protrusion (measured in μm) (Left) and on the number of protrusions (Right). Normalized histograms associated to in silico experiments (Bottom) based on the length of the longest protrusion (measured in μm) (Left) and on the number of protrusions (Right). In silico experiments were generated using one of the best parametrizations suggested by SigOpt with metrics BC_llp = 0.87 and BC_np = 0.81.

The values of both metrics BC_llp and BC_np for every suggested parametrization by SigOpt are shown in Figure 8. SigOpt is able to find parametrizations with higher values of the BC_llp (even higher than 0.9) than the BC_np (always lower than 0.8). Besides, the majority of the parametrizations suggested by SigOpt are higher than 0.7 for both metrics (35.67%), with slightly better results for the metric related with the length of the longest protrusion (Figure 8).

FIGURE 8

Figure 8. Values associated to both metrics (BC_llp and BC_np) for the 300 model parametrizations suggested by SigOpt during the calibration process. Red circles are associated to every parametrization tested whereas the blue ones represent Pareto optimal points (parametrizations where one metric cannot be improved without another metric suffering) and form an approximate Pareto frontier.

A total of nine parameters of the model have been calibrated (Table 5). The range of values for each parameter has also been established in order to define the search space (Table 5). Note that in this case, the search space includes both continuous regions in the real space and discrete values for integer parameters. Thus becoming a mixed-integer programming problem, much harder to be optimized than just real spaces (non-linear optimization) or integer spaces (combinatorial optimization). For some parameters we have established a range based on the values used in Ribeiro et al. (2017), whereas for others such as E_protrusion we have determined a range based on values found in literature. In addition, for the parameters related to s signal (s_birth, s_exp, s_ret, and s_binary), we have analyzed the values of s at different time steps. These ranges should be biologically relevant. For example, the range of the parameter E_protrusion (protrusions elastic modulus) includes the value given in Mofrad and Kamm (2006) and Li et al. (2014). We have also automatically discarded any parametrization with s_ret ≥ s_birth or s_ret ≥ s_exp or s_exp ≥ s_birth because from a biological perspective they are invalid (the minimal amount of signal required for the onset of new protrusions, s_birth, and for the reinforcement of pre-existing protrusions, s_exp, cannot be lower than the minimal amount of chemotactic signal s required to remain active and not disappear, s_ret; the minimal amount of signal required for the onset of new protrusions, s_birth, cannot be lower than the minimal amount for the reinforcement of pre-existing protrusions either). The parametrization selected as the optimal one after 300 iterations of the calibration process using SigOpt is summarized in Table 2. For example, the best value for the elastic modulus is 10⁷ Pa. The best parametrization, with metrics BC_llp = 0.87 and BC_np = 0.81, have been selected due to the balance between both metrics.

The advantage of having a probabilistic surrogate model of the metrics is that we can perform other types of data analysis during the optimization process. SigOpt also offers an importance analysis of each parameter on the metrics (see Figure 9), i.e., how influential each parameter is on the metrics, that is, how much the metric values change with variations of each parameter. This analysis gives us valuable insights on our model performance. Although every parameter has some influence over the metrics output, α_exp, a parameter which computes the free expansion/retraction stretch rate field (during the protrusion dynamics stage) is the most important parameter (24.06%). The parameters β_exp, also used to compute the free expansion/retraction stretch rate field, and s_binary, used to simplify the search of signal s peaks where protrusions centroids are localized, are the second and third most influential parameters on the metrics (15.25 and 13.76%, respectively). Lastly, s_birth and s_exp are the ones with the least importance on our evaluation metrics (4.76 and 3.86%, respectively).

FIGURE 9

Figure 9. Parameters sensitivity based on SigOpt analysis of each parameter importance on the proposed evaluation metrics.

Finally, we validate this computational model using qualitative observations based on cell motion of migrating cells surrounded by different chemoattractant gradients. Figure 10 (left) shows that as the PDGF gradient grows, cell's velocity in the direction of the chemotactic gradient increases too. Thus, cells are following a more directional trajectory which agrees with experimental observations from Bosgraaf and Van Haastert (2009). On the other hand, Figure 10 (right) shows that as the PDGF concentration surrounding the cell increases, cell's velocity in the direction of the chemotactic gradient decreases. In this case, cells are following a more random trajectory. This fall in the effective speed of the cell is thought to be associated with receptor saturation (Song et al., 2006).

FIGURE 10

Figure 10. Cell migration speed statistical analysis for 20 simulations using the parametrization selected during the calibration process and associated with six different extracellular environments. (Left) three of these extracellular environments include different PDGF gradients (10⁻¹, 10⁰, 10¹ μM/mm) but a fixed PDGF concentration at the initial cell's position (0.8 μM). (Right) three extracellular environments include a fixed PDGF gradient (10⁰ μM/mm) but different PDGF concentrations at the initial cell's position (0.08, 0.8, and 8.0 μM).

Discussion

Understanding the process of cell migration is a really difficult endeavor since it is a biological process coordinated by multiple factors. Temperature (Higazi et al., 1996), adhesion sites in the ECM (Cukierman et al., 2001), ECM mechanical properties and architecture (Wolf et al., 2013) as well as the gradient of chemical factors (Devreotes and Janetopoulos, 2003), modulate cell migration, by regulating the signaling pathways and intracellular cytoskeleton and adhesion organization (Paul et al., 2016).

According to our experimental observations (Moreno-Arotzena et al., 2015; Del Amo et al., 2017; Movilla et al., 2018) cells tend to present two different behaviors: they increase the number of stable protrusions, in which case each protrusion is shorter; or they decrease the number of stable protrusions, in which case at least some of them are longer. In the first case, protrusions compete and there is not any preferential movement. In the second case, normally cells present a defined movement in the direction of the longest protrusion.

Several assumptions are made regarding the mechanical model of the ECM. First, we consider the ECM as an isotropic material. Nevertheless, the ECM is anisotropic due to the different fiber directions (Valero et al., 2018). Second, the mechanical properties of the ECM are assumed homogeneous, thus we do not consider the heterogeneity associated to the distribution of the fibers. Third, ECM remodeling is not considered in this model. However, this is an acceptable approximation for preliminary studies of cell motility in collagen gels, which allows us to use the Eshelby's theory.

Due to the complexity of cell migration, computational models have been widely used to improve its understanding (Rangarajan and Zaman, 2008; Mak et al., 2016; Chen et al., 2017). Cell migration include several stochastic processes such as the evolution of chemically reacting systems. The Stochastic Simulation Algorithm (SSA) (Gillespie, 1976, 1977) has been widely used to numerically simulate the stochastic behavior of biochemical reactions. However, the SSA is considered too slow for many practical applications (Gillespie, 2001, 2007). This effect occurs clearly in our case: even though the SSA offers an exact solution, simulations take too long to finish (an average of 10.77 h of execution time for each simulations of 4 h of cell migration). The tau leaping algorithm has been considered a good fit for our purposes: it gives us a “good-enough” approximation (see footnote 1) of the temporal evolution of our biochemical system and allows us to optimize the numerical performance of our mechanochemical model (an average of 1.28 h of execution time for each simulation of 4 h of cell migration). Thus, reducing the computational cost to almost an order of magnitude.

Although in most computational works (Bauer et al., 2009; Bentley et al., 2009; Vermolen and Javierre, 2012; Daub and Merks, 2013; Talkenberger et al., 2017; Escribano et al., 2018; González-Valverde and García-Aznar, 2018; Kim et al., 2018; Moure and Gomez, 2018) authors perform strong efforts to validate models comparing experimental results with numerical ones, there is a lack of full integration of both kind of results. However, this paper presents a relevant step forward in this direction, showing a novel methodology that integrates both modeling strategies (in vitro and in silico) by means of the application of Bayesian optimization during the calibration process.

The complexity of the calibration process of any model grows rapidly with the number of parameters. Another factor that greatly increases the complexity of the calibration process is the stochastic nature of some biological models such as the one presented in this paper. Stochastic models require the execution of several simulations for each model parametrization in order to capture the results variation associated to the stochastic randomness. Moreover, if the execution of each simulation takes more than a couple of minutes, a manual approach for this calibration process becomes highly prone to inefficiencies.

When choosing the values for each model parameter using such a manual approach, it is usually the case that researchers turn to literature as their starting point. Then, they perform some manual tuning so simulations results fit approximately the experimental data. Generally, researchers start by modifying just a couple of parameters using some values considered biologically relevant. Then, they analyze how those parameters influence the model output based on the different values tested. They iterate over this process by picking a couple of the remaining parameters in every iteration—ideally, the selected parameters in each iteration are related to each other. This manual approach is really tedious since the modification of some parameters can potentially require the recalibration of some already calibrated parameters. If the model includes a large number of parameters, researchers could start this tuning process by performing a sensitivity analysis (Saltelli, 2002; Bauer et al., 2009; Bentley et al., 2009; Borau et al., 2012; Vermolen and Javierre, 2012; Daub and Merks, 2013; Escribano et al., 2015, 2018; Talkenberger et al., 2017) in order to focus on those parameters with a higher importance on the model output. Due to computational and time restrictions, this manual step does not generally include more than a couple of iterations, even though it is becoming more and more common to have access to a High-Throughput Computing (HTC) environment—which can reduce the required times to run those simulations by parallelizing them.

This paper proposes the application of the Bayesian optimization technique to reduce these inefficiencies. Bayesian optimization, which has been applied to solve a wide range of problems such as machine learning applications (Snoek et al., 2012), robot planning (Martinez-Cantin et al., 2009), simulation design (Brochu et al., 2010), biochemistry (Czarnecki et al., 2015), and dynamical modeling of biological systems (Ulmasov et al., 2016), offers an automated approach for this calibration process. Furthermore, the Bayesian optimization technique is able to minimize the number of parametrizations to test on the computational model and find a good enough fit to in vitro observations. In our case, from the 300 different parametrizations tested during the calibration process, only 6 parametrizations (2%) have the two metrics considered (BC_llp and BC_np) below 0.5. On the other hand, SigOpt suggests 107 parametrizations (35.67%) with both metrics above 0.7.

Clearly, the methodology here presented—based on the application of Bayesian optimization to compare the results of in vitro and in silico experiments—has allowed to identify the key parameters that regulates individual 3D fibroblast migration embedded in a collagen-based matrix. In particular, this novel methodology has been applied during the development of a stochastic model that simulates a chemically reacting system based on the biochemical interaction between the PDGF and a specific type of cell surface receptors, the RTKs. This interaction, in turn, triggers a metabolic cascade of internal signaling that activates a cellular chemosensing mechanism. Moreover, the model's calibration has been proven to be a valid and not an overfitted one during the final validation process. In order to validate the selected parametrization, we have simulated cell migration with a diffused chemoattractant factor throughout the ECM and qualitatively compare observations based on cell's velocity in the direction of the chemotactic gradient with results from previous experimental works (Song et al., 2006; Bosgraaf and Van Haastert, 2009). Our results are in agreement with those from in vitro experiments, cells follow a more directional motion as the chemoattractant gradient increases. However, when the chemoattractant concentration surrounding the cell reaches a saturation point cells start to lose the ability to sense the chemical cues.

In conclusion, the tau leaping algorithm allows to optimize the performance of stochastic models based on biochemical kinetics by greatly reducing the execution time of its simulations. In addition, by means of Bayesian optimization it is possible to perform model parameters calibration in a very efficient and completely automatic way. As a result, this novel methodology will improve the development of in silico models for a better understanding of cell migration.

Author Contributions

FM-C, MG-B, and JG-A designed research; FM-C, MG-B, YJ-L, and JG-A performed research; YJ-L performed in vitro experiments, FM-C, MG-B, RM-C, and JG-A analyzed data; FM-C, MG-B, YJ-L, RM-C, and JG-A wrote the paper; FM-C, MG-B, and RM-C defined Bayesian optimization setup and FM-C built a code for the model and performed all simulations.

Funding

FM-C was supported by Spanish Ministry of Economy and Competitiveness (Grant no: BES-2016-076291). MG-B, YJ-L, and JG-A were supported by the European Research Council (Grant no: ERC2012-StG306571) and the Spanish Ministry of Economy and Competitiveness (Grant no: DPI2015-64221-C2-1-R). RM-C was supported by the Spanish Ministry of Economy and Competitiveness (Grant no: DPI2015-65962-R).

Conflict of Interest Statement

RM-C was employed by company SigOpt, Inc.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer AL and handling editor declared their shared affiliation at the time of the review.

Acknowledgments

The authors would like to acknowledge Frederico O. Ribeiro assistance during the model's development and SigOpt Inc. for their assistance and support with their Bayesian optimization framework.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2018.01246/full#supplementary-material

Footnotes

1. ^The “good-enough” expression used here to describe the accuracy of the tau-leaping algorithm comes from previous works such as Lok (2004) where he states that “One acceleration strategy is to abandon absolute mathematical precision in favor of a good-enough approximation. Gillespie has also been a pioneer in this effort. One of his strategies is called ‘tau-leaping'.” This statement is considered valid as long as the leap condition is satisfied, i.e., as long as the probability of each reaction taking place does not change significantly over the time leap.

2. ^The leap condition is an accuracy-assuring restriction which states that during the time interval [t, t+τ] the probability of each reaction channel R_j being fired should remain approximately constant even though all reaction channels may be fired several times.

3. ^https://sigopt.com/

References

Alber, M., Chen, N., Lushnikov, P. M., and Newman, S. A. (2007). continuous macroscopic limit of a discrete stochastic model for interaction of living cells. Phys. Rev. Lett. 99:168102. doi: 10.1103/PhysRevLett.99.168102

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, D. F., Ganguly, A., and Kurtz, T. G. (2011). Error analysis of tau-leap simulation methods. Ann. Appl. Probab. 21, 2226–2262. doi: 10.1214/10-AAP756

CrossRef Full Text | Google Scholar

Bauer, A. L., Jackson, T. L., and Jiang, Y. (2009). Topography of extracellular matrix mediates vascular morphogenesis and migration speeds in angiogenesis. PLoS Comput. Biol. 5:e1000445. doi: 10.1371/journal.pcbi.1000445

PubMed Abstract | CrossRef Full Text | Google Scholar

Bazmara, H., Soltani, M., Sefidgar, M., Bazargan, M., Mousavi Naeenian, M., and Rahmim, A. (2015). The vital role of blood flow-induced proliferation and migration in capillary network formation in a multiscale model of angiogenesis. PLoS ONE 10:e0128878. doi: 10.1371/journal.pone.0128878

PubMed Abstract | CrossRef Full Text | Google Scholar

Bear, J. E., and Haugh, J. M. (2014). Directed migration of mesenchymal cells: where signaling and the cytoskeleton meet. Curr. Opin. Cell Biol. 30, 74–82. doi: 10.1016/J.CEB.2014.06.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Beauchemin, C., Dixit, N. M., and Perelson, A. S. (2007). Characterizing T cell movement within lymph nodes in the absence of antigen. J. Immunol. 178, 5505–5512. doi: 10.4049/JIMMUNOL.178.9.5505

PubMed Abstract | CrossRef Full Text | Google Scholar

Bentley, K., Mariggi, G., Gerhardt, H., and Bates, P. A. (2009). Tipping the balance: robustness of tip cell selection, migration and fusion in angiogenesis. PLoS Comput. Biol. 5:e1000549. doi: 10.1371/journal.pcbi.1000549

PubMed Abstract | CrossRef Full Text | Google Scholar

Bogle, G., and Dunbar, P. R. (2010). T cell responses in lymph nodes. Wiley Interdiscip. Rev. Syst. Biol. Med. 2, 107–116. doi: 10.1002/wsbm.47

PubMed Abstract | CrossRef Full Text | Google Scholar

Borau, C., Kamm, R. D., and García-Aznar, J. M. (2011). Mechano-sensing and cell migration: a 3D model approach. Phys. Biol. Phys. Biol. Phys. Biol 8, 66008–66013. doi: 10.1088/1478-3975/8/6/066008

PubMed Abstract | CrossRef Full Text | Google Scholar

Borau, C., Kim, T., Bidone, T., García-Aznar, J. M., and Kamm, R. D. (2012). Dynamic mechanisms of cell rigidity sensing: insights from a computational model of actomyosin networks. PLoS ONE 7:e49174. doi: 10.1371/journal.pone.0049174

PubMed Abstract | CrossRef Full Text | Google Scholar

Bosgraaf, L., and Van Haastert, P. J. M. (2009). Navigation of chemotactic cells by parallel signaling to pseudopod persistence and orientation. PLoS ONE 4:e6842. doi: 10.1371/journal.pone.0006842

PubMed Abstract | CrossRef Full Text | Google Scholar

Brochu, E., Brochu, T., and de Freitas, N. (2010). “A Bayesian interactive optimization approach to procedural animation design,” in Proc. 2010 ACM SIGGRAPH/Eurographics Symp. Comput. Animat. (Madrid), 103–112.

Google Scholar

Camley, B. A., and Rappel, W.-J. (2017). Physical models of collective cell motility: from cell to tissue. J. Phys. D Appl. Phys. 50:113002. doi: 10.1088/1361-6463/aa56fe

PubMed Abstract | CrossRef Full Text | Google Scholar

Campellone, K. G., and Welch, M. D. (2010). A nucleator arms race: cellular control of actin assembly. Nat. Rev. Mol. Cell Biol. 11, 237–251. doi: 10.1038/nrm2867

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, R., Björndahl, M. A., Religa, P., Clasper, S., Garvin, S., Galter, D., et al. (2004). PDGF-BB induces intratumoral lymphangiogenesis and promotes lymphatic metastasis. Cancer Cell 6, 333–345. doi: 10.1016/J.CCR.2004.08.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, Y., Gillespie, D. T., and Petzold, L. R. (2006). Efficient step size selection for the tau-leaping simulation method. J. Chem. Phys. 124:044109. doi: 10.1063/1.2159468

PubMed Abstract | CrossRef Full Text | Google Scholar

Cazzaniga, P., Pescini, D., Besozzi, D., and Mauri, G. (2006). Tau Leaping Stochastic Simulation Method in P Systems. Berlin; Heidelberg: Springer Berlin Heidelberg.

Google Scholar

Chen, J., Weihs, D., and Vermolen, F. J. (2017). A model for cell migration in non-isotropic fibrin networks with an application to pancreatic tumor islets. Biomech. Model. Mechanobiol. 17, 367–386. doi: 10.1007/s10237-017-0966-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, R. R., Silva, E. A., Yuen, W. W., and Mooney, D. J. (2007). Spatio–temporal VEGF and PDGF delivery patterns blood vessel formation and maturation. Pharm. Res. 24, 258–264. doi: 10.1007/s11095-006-9173-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Comaniciu, D., Ramesh, V., and Meer, P. (2000). “Real-time tracking of non-rigid objects using mean shift,” in Proceedings IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2000 (Cat. No.PR00662) (Hilton Head Island, SC: IEEE Comput. Soc), 142–149.

Google Scholar

Condeelis, J., and Pollard, J. W. (2006). Macrophages: obligate partners for tumor cell migration, invasion, and metastasis. Cell 124, 263–266. doi: 10.1016/J.CELL.2006.01.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Condeelis, J., Singer, R. H., and Segall, J. E. (2005). THE GREAT ESCAPE: when cancer cells hijack the genes for chemotaxis and motility. Annu. Rev. Cell Dev. Biol. 21, 695–718. doi: 10.1146/annurev.cellbio.21.122303.120306

PubMed Abstract | CrossRef Full Text | Google Scholar

Cukierman, E., Pankov, R., Stevens, D. R., and Yamada, K. M. (2001). Taking cell-matrix adhesions to the third dimension. Science 294, 1708–1712. doi: 10.1126/science.1064829

PubMed Abstract | CrossRef Full Text | Google Scholar

Czarnecki, W. M., Podlewska, S., and Bojarski, A. J. (2015). Robust optimization of SVM hyperparameters in the classification of bioactive compounds. J. Cheminform. 7:38. doi: 10.1186/s13321-015-0088-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Daub, J. T., and Merks, R. M. H. (2013). A cell-based model of extracellular-matrix-guided endothelial cell migration during angiogenesis. Bull. Math. Biol. 75, 1377–1399. doi: 10.1007/s11538-013-9826-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Del Amo, C., Borau, C., Movilla, N., Asín, J., and García-Aznar, J. M. (2017). Quantifying 3D chemotaxis in microfluidic-based chips with step gradients of collagen hydrogel concentrations. Integr. Biol. 9, 339–349. doi: 10.1039/C7IB00022G

PubMed Abstract | CrossRef Full Text | Google Scholar

Devreotes, P., and Janetopoulos, C. (2003). Eukaryotic chemotaxis: distinctions between directional sensing and polarization. J. Biol. Chem. 278, 20445–20448. doi: 10.1074/jbc.R300010200

PubMed Abstract | CrossRef Full Text | Google Scholar

Elangovan, S., D'Mello, S. R., Hong, L., Ross, R. D., Allamargot, C., Dawson, D. V., et al. (2014). The enhancement of bone regeneration by gene activated matrix encoding for platelet derived growth factor. Biomaterials 35, 737–747. doi: 10.1016/J.BIOMATERIALS.2013.10.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Escribano, J., Sánchez, M. T., and García-Aznar, J. M. (2015). Modeling the formation of cell-matrix adhesions on a single 3D matrix fiber. J. Theor. Biol. 384, 84–94. doi: 10.1016/J.JTBI.2015.07.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Escribano, J., Sunyer, R., Sánchez, M. T., Trepat, X., Roca-Cusachs, P., and García-Aznar, J. M. (2018). A hybrid computational model for collective cell durotaxis. Biomech. Model. Mechanobiol. 17, 1037–1052. doi: 10.1007/s10237-018-1010-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Eshelby, J. D. (1957). The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc. R Soc. A Math. Phys. Eng. Sci. 241, 376–396. doi: 10.1098/rspa.1957.0133

CrossRef Full Text | Google Scholar

Fraley, S., Feng, Y., Krishnamurthy, R., Kim, D. H., Celedon, A., Longmore, G. D. D., et al. (2010). A distinctive role for focal adhesion proteins in three-dimensional cell motility. Nat. Cell Biol. 12, 598–604. doi: 10.1038/ncb2062

PubMed Abstract | CrossRef Full Text | Google Scholar

Franz, C. M., Jones, G. E., and Ridley, A. J. (2002). Cell migration in development and disease. Dev. Cell 2, 153–158. doi: 10.1016/S1534-5807(02)00120-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedl, P., and Wolf, K. (2010). Plasticity of cell migration: a multiscale tuning model. J. Cell Biol. 188, 11–19. doi: 10.1083/jcb.200909003

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedlaender, G. E., Lin, S., Solchaga, L. A., Snel, L. B., and Lynch, S. E. (2013). The role of recombinant human platelet-derived growth factor-BB (rhPDGF-BB) in orthopaedic bone repair and regeneration. Curr. Pharm. Des. 19, 3384–3390. doi: 10.2174/1381612811319190005

PubMed Abstract | CrossRef Full Text | Google Scholar

Gillespie, D. T. (1976). A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys. 22, 403–434. doi: 10.1016/0021-9991(76)90041-3

CrossRef Full Text | Google Scholar

Gillespie, D. T. (1977). Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 81, 2340–2361. doi: 10.1021/j100540a008

CrossRef Full Text | Google Scholar

Gillespie, D. T. (2001). Approximate accelerated stochastic simulation of chemically reacting systems. J. Chem. Phys. 115, 1716–1733. doi: 10.1063/1.1378322

CrossRef Full Text | Google Scholar

Gillespie, D. T. (2007). Stochastic simulation of chemical kinetics. Annu. Rev. Phys. Chem. 58, 35–55. doi: 10.1146/annurev.physchem.58.032806.104637

PubMed Abstract | CrossRef Full Text | Google Scholar

González-Valverde, I., and García-Aznar, J. M. (2018). Mechanical modeling of collective cell migration: an agent-based and continuum material approach. Comput. Methods Appl. Mech. Eng. 337, 246–262. doi: 10.1016/J.CMA.2018.03.036

CrossRef Full Text | Google Scholar

González-Valverde, I., Semino, C., and García-Aznar, J. M. (2016). Phenomenological modelling and simulation of cell clusters in 3D cultures. Comput. Biol. Med. 77, 249–260. doi: 10.1016/J.COMPBIOMED.2016.08.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansen, N., Müller, S. D., and Koumoutsakos, P. (2003). Reducing the time complexity of the derandomized evolution strategy with Covariance Matrix Adaptation (CMA-ES). Evol. Comput. 11, 1–18. doi: 10.1162/106365603321828970

PubMed Abstract | CrossRef Full Text | Google Scholar

Hatakeyama, M., Kimura, S., Naka, T., Kawasaki, T., Yumoto, N., Ichikawa, M., et al. (2003). A computational model on the modulation of mitogen-activated protein kinase (MAPK) and Akt pathways in heregulin-induced ErbB signalling. Biochem. J. 373, 451–463. doi: 10.1042/BJ20021824

PubMed Abstract | CrossRef Full Text | Google Scholar

Hawkins, P. T., Anderson, K. E., Davidson, K., and Stephens, L. R. (2006). Signalling through Class I PI3Ks in mammalian cells. Biochem. Soc. Trans. 34, 647–662. doi: 10.1042/BST0340647

PubMed Abstract | CrossRef Full Text | Google Scholar

Heinecke, K., Seher, A., Schmitz, W., Mueller, T. D., Sebald, W., and Nickel, J. (2009). Receptor oligomerization and beyond: a case study in bone morphogenetic proteins. BMC Biol. 7:59. doi: 10.1186/1741-7007-7-59

PubMed Abstract | CrossRef Full Text | Google Scholar

Higazi, A. A., Kniss, D., Manuppello, J., Barnathan, E. S., and Cines, D. B. (1996). Thermotaxis of human trophoblastic cells. Placenta 17, 683–687. doi: 10.1016/S0143-4004(96)80019-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Jilkine, A., and Edelstein-Keshet, L. (2011). A Comparison of mathematical models for polarization of single eukaryotic cells in response to guided cues. PLoS Comput. Biol. 7:e1001121. doi: 10.1371/journal.pcbi.1001121

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, D. R., Schonlau, M., and Welch, W. J. (1998). Efficient global optimization of expensive black-box functions. J. Glob. Optim. 13, 455–492. doi: 10.1023/A:1008306431147

CrossRef Full Text | Google Scholar

Jones, E., Oliphant, T., and Peterson, P. (2001). {SciPy}: Open Source Scientific Tools for {Python}. Available online at: https://www.scipy.org/citing.html

Kass, R. E., and Raftery, A. E. (1995). Bayes factors. J. Am. Stat. Assoc. 90:773. doi: 10.2307/2291091

CrossRef Full Text | Google Scholar

Kennedy, J., and Eberhart, R. (1995). “Particle swarm optimization,” in Proceedings of IEEE International Conference on Neural Networks (ICNN) (Nagoya: IEEE), 1942–1948.

Kim, M. C., Silberberg, Y. R., Abeyaratne, R., Kamm, R. D., and Asada, H. H. (2018). Computational modeling of three-dimensional ECM-rigidity sensing to guide directed cell migration. Proc. Natl. Acad. Sci. U.S.A. 115, E390–E399. doi: 10.1073/pnas.1717230115

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, M. C., Whisler, J., Silberberg, Y. R., Kamm, R. D., and Asada, H. H. (2015). Cell invasion dynamics into a three dimensional extracellular matrix fibre network. PLOS Comput. Biol. 11:e1004535. doi: 10.1371/journal.pcbi.1004535

PubMed Abstract | CrossRef Full Text | Google Scholar

Knecht, A. K., and Bronner-Fraser, M. (2002). Induction of the neural crest: a multigene process. Nat. Rev. Genet. 3, 453–461. doi: 10.1038/nrg819

PubMed Abstract | CrossRef Full Text | Google Scholar

Knowles, J. (2006). ParEGO: a hybrid algorithm with on-line landscape approximation for expensive multiobjective optimization problems. IEEE Trans. Evol. Comput. 10, 50–66. doi: 10.1109/TEVC.2005.851274

CrossRef Full Text | Google Scholar

Lamalice, L., Le Boeuf, F., and Huot, J. (2007). Endothelial cell migration during angiogenesis. Circ. Res. 100, 782–794. doi: 10.1161/01.RES.0000259593.07661.1e

PubMed Abstract | CrossRef Full Text | Google Scholar

Lämmermann, T., Bader, B. L., Monkley, S. J., Worbs, T., Wedlich-Söldner, R., Hirsch, K., et al. (2008). Rapid leukocyte migration by integrin-independent flowing and squeezing. Nature 453, 51–55. doi: 10.1038/nature06887

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, T., Gu, Y. T., Oloyede, A., and Yarlagadda, P. K. (2014). Molecular investigation of the mechanical properties of single actin filaments based on vibration analyses. Comput. Methods Biomech. Biomed. Eng. 17, 616–622. doi: 10.1080/10255842.2012.706279

PubMed Abstract | CrossRef Full Text | Google Scholar

Liou, Y. R., Torng, W., Kao, Y. C., Sung, K. B., Lee, C. H., and Kuo, P. L. (2014). Substrate stiffness regulates filopodial activities in lung cancer cells. PLoS ONE 9:e89767. doi: 10.1371/journal.pone.0089767

PubMed Abstract | CrossRef Full Text | Google Scholar

Lok, L. (2004). The need for speed in stochastic simulation. Nat. Biotechnol. 22, 964–965. doi: 10.1038/nbt0804-964

PubMed Abstract | CrossRef Full Text | Google Scholar

Luster, A. D., Alon, R., and von Andrian, U. H. (2005). Immune cell migration in inflammation: present and future therapeutic targets. Nat. Immunol. 6, 1182–1190. doi: 10.1038/ni1275

PubMed Abstract | CrossRef Full Text | Google Scholar

Mak, M., Spill, F., Kamm, R. D., and Zaman, M. H. (2016). Single-cell migration in complex microenvironments: mechanics and signaling dynamics. J. Biomech. Eng. 138:021004. doi: 10.1115/1.4032188

PubMed Abstract | CrossRef Full Text | Google Scholar

Mark, C., Metzner, C., Lautscham, L., Strissel, P. L., Strick, R., and Fabry, B. (2018). Bayesian model selection for complex dynamic systems. Nat. Commun. 9:1803. doi: 10.1038/s41467-018-04241-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, P., and Parkhurst, S. M. (2004). Parallels between tissue repair and embryo morphogenesis. Development 131, 3021–3034. doi: 10.1242/dev.01253

PubMed Abstract | CrossRef Full Text | Google Scholar

Martinez-Cantin, R. (2014). BayesOpt: a bayesian optimization library for nonlinear optimization, experimental design and bandits. J. Mach. Learn. Res. 15, 3915–3919. Available online at: http://www.jmlr.org/papers/v15/martinezcantin14a.html

Google Scholar

Martinez-Cantin, R., de Freitas, N., Brochu, E., Castellanos, J., and Doucet, A. (2009). A Bayesian exploration-exploitation approach for optimal online sensing and planning with a visually guided mobile robot. Auton. Robots 27, 93–103. doi: 10.1007/s10514-009-9130-2

CrossRef Full Text | Google Scholar

Martinez-Cantin, R., Tee, K., and McCourt, M. (2018). “Practical Bayesian optimization in the presence of outliers,” in International Conference on Artificial Intelligence and Statistics (AISTATS) (Lanzarote).

Google Scholar

Meineke, F. A., Potten, C. S., and Loeffler, M. (2001). Cell migration and organization in the intestinal crypt using a lattice-free model. Cell Prolif. 34, 253–266. doi: 10.1046/j.0960-7722.2001.00216.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Milde, F., Tauriello, G., Haberkern, H., and Koumoutsakos, P. (2014). SEM++: a particle model of cellular growth, signaling and migration. Comput. Particle Mech. 1, 211–227. doi: 10.1007/s40571-014-0017-4

CrossRef Full Text | Google Scholar

Mofrad, M., and Kamm, R. (eds.). (2006). Cytoskeletal Mechanics: Models and Measurements in Cell Mechanics (Cambridge Texts in Biomedical Engineering). Cambridge: Cambridge University Press.

Moreno-Arotzena, O., Borau, C., Movilla, N., Vicente-Manzanares, M., and García-Aznar, J. M. (2015). Fibroblast migration in 3D is controlled by haptotaxis in a non-muscle myosin II-dependent manner. Ann. Biomed. Eng. 43, 3025–3039. doi: 10.1007/s10439-015-1343-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Moreno-Arotzena, O., Mendoza, G., Cóndor, M., Rüberg, T., and García-Aznar, J. M. (2014). Inducing chemotactic and haptotactic cues in microfluidic devices for three-dimensional in vitro assays. Biomicrofluidics 8:064122. doi: 10.1063/1.4903948

PubMed Abstract | CrossRef Full Text | Google Scholar

Moure, A., and Gomez, H. (2017). Phase-field model of cellular migration: three-dimensional simulations in fibrous networks. Comput. Methods Appl. Mech. Eng. 320, 162–197. doi: 10.1016/J.CMA.2017.03.025

CrossRef Full Text | Google Scholar

Moure, A., and Gomez, H. (2018). Three-dimensional simulation of obstacle-mediated chemotaxis. Biomech. Model. Mechanobiol. doi: 10.1007/s10237-018-1023-x. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

Movilla, N., Borau, C., Valero, C., and García-Aznar, J. M. (2018). Degradation of extracellular matrix regulates osteoblast migration: a microfluidic-based study. Bone 107, 10–17. doi: 10.1016/j.bone.2017.10.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Nielsen, F. (2014). Generalized Bhattacharyya and Chernoff upper bounds on Bayes error using quasi-arithmetic means. Pattern Recognit. Lett. 42, 25–34. doi: 10.1016/J.PATREC.2014.01.002

CrossRef Full Text | Google Scholar

Nocedal, J. (1980). Updating quasi-Newton matrices with limited storage. Math. Comput. 35, 773–773. doi: 10.1090/S0025-5718-1980-0572855-7

CrossRef Full Text | Google Scholar

Norton, K. A., and Popel, A. S. (2016). Effects of endothelial cell proliferation and migration rates in a computational model of sprouting angiogenesis. Sci. Rep. 6:36992. doi: 10.1038/srep36992

PubMed Abstract | CrossRef Full Text | Google Scholar

Paralkar, V. M., Weeks, B. S., Yu, Y. M., Kleinman, H. K., and Reddi, A. H. (1992). Recombinant human bone morphogenetic protein 2B stimulates PC12 cell differentiation: potentiation and binding to type IV collagen. J. Cell Biol. 119, 1721–1728. doi: 10.1083/JCB.119.6.1721

PubMed Abstract | CrossRef Full Text | Google Scholar

Paul, C. D., Hung, W. C., Wirtz, D., and Konstantopoulos, K. (2016). Engineered models of confined cell migration. Annu. Rev. Biomed. Eng. 18, 159–180. doi: 10.1146/annurev-bioeng-071114-040654

PubMed Abstract | CrossRef Full Text | Google Scholar

Poukkula, M., Cliffe, A., Changede, R., and Rørth, P. (2011). Cell behaviors regulated by guidance cues in collective migration of border cells. J. Cell Biol. 192, 513–524. doi: 10.1083/jcb.201010003

PubMed Abstract | CrossRef Full Text | Google Scholar

Provenzano, P. P., Inman, D. R., Eliceiri, K. W., Trier, S. M., and Keely, P. J. (2008). Contact guidance mediated three-dimensional cell migration is regulated by Rho/ROCK-dependent matrix reorganization. Biophys. J. 95, 5374–5384. doi: 10.1529/BIOPHYSJ.108.133116

PubMed Abstract | CrossRef Full Text | Google Scholar

Rangarajan, R., and Zaman, M. H. (2008). Modeling cell migration in 3D. Cell Adh. Migr. 2, 106–109. doi: 10.4161/cam.2.2.6211

PubMed Abstract | CrossRef Full Text | Google Scholar

Reina-Romo, E., Gómez-Benito, M. J., Domínguez, J., and García-Aznar, J. M. (2012). A lattice-based approach to model distraction osteogenesis. J. Biomech. 45, 2736–2742. doi: 10.1016/j.jbiomech.2012.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Ribeiro, F. O., Gómez-Benito, M. J., Folgado, J., Fernandes, P. R., and García-Aznar, J. M. (2017). Computational model of mesenchymal migration in 3D under chemotaxis. Comput. Methods Biomech. Biomed. Eng. 20, 59–74. doi: 10.1080/10255842.2016.1198784

PubMed Abstract | CrossRef Full Text | Google Scholar

Roca-Cusachs, P., Sunyer, R., and Trepat, X. (2013). Mechanical guidance of cell migration: lessons from chemotaxis. Curr. Opin. Cell Biol. 25, 543–549. doi: 10.1016/J.CEB.2013.04.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Saltelli, A. (2002). Sensitivity analysis for importance assessment. Risk Anal. 22, 579–590. doi: 10.1111/0272-4332.00040

PubMed Abstract | CrossRef Full Text | Google Scholar

Saltzman, E. A., Drew, J. H., Leemis, L. M., and Henderson, S. G. (2012). Simulating multivariate nonhomogeneous poisson processes using projections. ACM Trans. Model. Comput. Simul. 22, 1–13. doi: 10.1145/2331140.2331143

CrossRef Full Text | Google Scholar

Schlüter, D. K., Ramis-Conde, I., and Chaplain, M. A. (2012). Computational modeling of single-cell migration: the leading role of extracellular matrix fibers. Biophys. J. 103, 1141–1151. doi: 10.1016/J.BPJ.2012.07.048

PubMed Abstract | CrossRef Full Text | Google Scholar

Scianna, M., and Preziosi, L. (2014). A cellular Potts model for the MMP-dependent and -independent cancer cell migration in matrix microtracks of different dimensions. Comput. Mech. 53, 485–497. doi: 10.1007/s00466-013-0944-6

CrossRef Full Text | Google Scholar

Scianna, M., Preziosi, L., and Wolf, K. (2012). A Cellular Potts model simulating cell migration on and in matrix environments. Math. Biosci. Eng. 10, 235–261. doi: 10.3934/mbe.2013.10.235

PubMed Abstract | CrossRef Full Text | Google Scholar

Serrano-Alcalde, F., García-Aznar, J. M., and Gómez-Benito, M. J. (2017). The role of nuclear mechanics in cell deformation under creeping flows. J. Theor. Biol. 432, 25–32. doi: 10.1016/J.JTBI.2017.07.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Shah, P., Keppler, L., and Rutkowski, J. (2014). A review of platelet derived growth factor playing pivotal role in bone regeneration. J. Oral Implantol. 40, 330–340. doi: 10.1563/AAID-JOI-D-11-00173

PubMed Abstract | CrossRef Full Text | Google Scholar

Shahriari, B., Swersky, K., Wang, Z., Adams, R. P., and de Freitas, N. (2016). Taking the human out of the loop: a review of bayesian optimization. Proc. IEEE 104, 148–175. doi: 10.1109/JPROC.2015.2494218

CrossRef Full Text | Google Scholar

Shaw, T. J., and Martin, P. (2009). Wound repair at a glance. J. Cell Sci. 122, 3209–3213. doi: 10.1242/jcs.031187

PubMed Abstract | CrossRef Full Text | Google Scholar

Snoek, J., Larochelle, H., and Adams, R. P. (2012). “Practical Bayesian optimization of machine learning algorithms,” in Conference on Neural Information Processing Systems (NIPS) (Lake Tahoe, NV).

Google Scholar

Song, L., Nadkarni, S. M., Bödeker, H. U., Beta, C., Bae, A., Franck, C., et al. (2006). Dictyostelium discoideum chemotaxis: threshold for directed motion. Eur. J. Cell Biol. 85, 981–989. doi: 10.1016/J.EJCB.2006.01.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Spill, F., Guerrero, P., Alarcon, T., Maini, P. K., and Byrne, H. M. (2015). Mesoscopic and continuum modelling of angiogenesis. J. Math. Biol. 70, 485–532. doi: 10.1007/s00285-014-0771-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Starke, J., Maaser, K., Wehrle-Haller, B., and Friedl, P. (2013). Mechanotransduction of mesenchymal melanoma cell invasion into 3D collagen lattices: filopod-mediated extension–relaxation cycles and force anisotropy. Exp. Cell Res. 319, 2424–2433. doi: 10.1016/J.YEXCR.2013.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, M., and Zaman, M. H. (2017). Modeling, signaling and cytoskeleton dynamics: integrated modeling-experimental frameworks in cell migration. Wiley Interdiscip. Rev. Syst. Biol. Med. 9:e1365. doi: 10.1002/wsbm.1365

PubMed Abstract | CrossRef Full Text | Google Scholar

Sunyer, R., Conte, V., Escribano, J., Elosegui-Artola, A., Labernadie, A., Valon, L., et al. (2016). Collective cell durotaxis emerges from long-range intercellular force transmission. Science 353, 1157–1161. doi: 10.1126/science.aaf7119

PubMed Abstract | CrossRef Full Text | Google Scholar

Swaney, K. F., Huang, C. H., and Devreotes, P. N. (2010). Eukaryotic chemotaxis: a network of signaling pathways controls motility, directional sensing, and polarity. Annu. Rev. Biophys. 39, 265–289. doi: 10.1146/annurev.biophys.093008.131228

PubMed Abstract | CrossRef Full Text | Google Scholar

Talkenberger, K., Cavalcanti-Adam, E. A., Voss-Böhme, A., and Deutsch, A. (2017). Amoeboid-mesenchymal migration plasticity promotes invasion only in complex heterogeneous microenvironments. Sci. Rep. 7:9237. doi: 10.1038/s41598-017-09300-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Te Boekhorst, V., Preziosi, L., and Friedl, P. (2016). Plasticity of cell migration in vivo and in silico. Annu. Rev. Cell Dev. Biol. 32, 491–526. doi: 10.1146/annurev-cellbio-111315-125201

PubMed Abstract | CrossRef Full Text | Google Scholar

Trichet, L., Le Digabel, J., Hawkins, R. J., Vedula, S. R., Gupta, M., Ribrault, C., et al. (2012). Evidence of a large-scale mechanosensing mechanism for cellular adaptation to substrate stiffness. Proc. Natl. Acad. Sci. U.S.A. 109, 6933–6938. doi: 10.1073/pnas.1117810109

PubMed Abstract | CrossRef Full Text | Google Scholar

Ueda, M., and Shibata, T. (2007). Stochastic signal processing and transduction in chemotactic response of eukaryotic cells. Biophys. J. 93, 11–20. doi: 10.1529/BIOPHYSJ.106.100263

PubMed Abstract | CrossRef Full Text | Google Scholar

Ulmasov, D., Baroukh, C., Chachuat, B., Deisenroth, M. P., and Misener, R. (2016). Bayesian optimization with dimension scheduling: application to biological systems. Comput. Aided Chem. Eng. 38, 1051–1056. doi: 10.1016/B978-0-444-63428-3.50180-6

CrossRef Full Text | Google Scholar

Valero, C., Amaveda, H., Mora, M., and García-Aznar, J. M. (2018). Combined experimental and computational characterization of crosslinked collagen-based hydrogels. PLoS ONE 13:e0195820. doi: 10.1371/journal.pone.0195820

PubMed Abstract | CrossRef Full Text | Google Scholar

Valero, C., Javierre, E., García-Aznar, J. M., and Gómez-Benito, M. J. (2014a). A cell-regulatory mechanism involving feedback between contraction and tissue formation guides wound healing progression. PLoS ONE 9:e92774. doi: 10.1371/journal.pone.0092774

PubMed Abstract | CrossRef Full Text | Google Scholar

Valero, C., Javierre, E., García-Aznar, J. M., and Gómez-Benito, M. J. (2014b). Nonlinear finite element simulations of injuries with free boundaries: application to surgical wounds. Int. J. Numer. Methods Biomed. Eng. 30, 616–633. doi: 10.1002/cnm.2621

PubMed Abstract | CrossRef Full Text | Google Scholar

van der Walt, S., Colbert, S. C., and Varoquaux, G. (2011). The NumPy array: a structure for efficient numerical computation. Comput. Sci. Eng. 13, 22–30. doi: 10.1109/MCSE.2011.37

CrossRef Full Text | Google Scholar

Van Liedekerke, P., Palm, M. M., Jagiella, N., and Drasdo, D. (2015). Simulating tissue mechanics with agent-based models: concepts, perspectives and some novel results. Comput. Particle Mech. 2, 401–444. doi: 10.1007/s40571-015-0082-3

CrossRef Full Text | Google Scholar

Vermolen, F. J., and Javierre, E. (2012). A finite-element model for healing of cutaneous wounds combining contraction, angiogenesis and closure. J. Math. Biol. 65, 967–996. doi: 10.1007/s00285-011-0487-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Weiger, M. C., Ahmed, S., Welf, E. S., and Haugh, J. M. (2010). Directional persistence of cell migration coincides with stability of asymmetric intracellular signaling. Biophys. J. 98, 67–75. doi: 10.1016/J.BPJ.2009.09.051

PubMed Abstract | CrossRef Full Text | Google Scholar

Wolf, K., Te Lindert, M., Krause, M., Alexander, S., Te Riet, J., Willis, A. L., et al. (2013). Physical limits of cell migration: control by ECM space and nuclear deformation and tuning by proteolysis and traction force. J. Cell Biol. 201, 1069–1084. doi: 10.1083/jcb.201210152

PubMed Abstract | CrossRef Full Text | Google Scholar

Zaman, M. H., Kamm, R. D., Matsudaira, P., and Lauffenburger, D. A. (2005). Computational model for cell migration in three-dimensional matrices. Biophys. J. 89, 1389–1397. doi: 10.1529/BIOPHYSJ.105.060723

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, J., and Mogilner, A. (2016). Comparison of cell migration mechanical strategies in three-dimensional matrices: a computational study. Interface Focus 6:20160040. doi: 10.1098/rsfs.2016.0040

PubMed Abstract | CrossRef Full Text | Google Scholar