Understanding glioblastoma invasion using physically-guided neural networks with internal variables.

Microfluidic capacities for both recreating and monitoring cell cultures have opened the door to the use of Data Science and Machine Learning tools for understanding and simulating tumor evolution under controlled conditions. In this work, we show how these techniques could be applied to study Glioblastoma, the deadliest and most frequent primary brain tumor. In particular, we study Glioblastoma invasion using the recent concept of Physically-Guided Neural Networks with Internal Variables (PGNNIV), able to combine data obtained from microfluidic devices and some physical knowledge governing the tumor evolution. The physics is introduced in the network structure by means of a nonlinear advection-diffusion-reaction partial differential equation that models the Glioblastoma evolution. On the other hand, multilayer perceptrons combined with a nodal deconvolution technique are used for learning the go or grow metabolic behavior which characterises the Glioblastoma invasion. The PGNNIV is here trained using synthetic data obtained from in silico tests created under different oxygenation conditions, using a previously validated model. The unravelling capacity of PGNNIV enables discovering complex metabolic processes in a non-parametric way, thus giving explanatory capacity to the networks, and, as a consequence, surpassing the predictive power of any parametric approach and for any kind of stimulus. Besides, the possibility of working, for a particular tumor, with different boundary and initial conditions, permits the use of PGNNIV for defining virtual therapies and for drug design, thus making the first steps towards in silico personalised medicine.

This is a PLOS Computational Biology Methods paper.

Introduction

Cancer is the second leading cause of death in the world, according to the World Health Organisation, and is responsible for about 10 million deaths per year. These figures are expected to rise up to 16 million deaths in 2040 [1]. Among the more than 200 types of cancer, Glioblastoma (GBM) is the most aggressive primary brain tumor, with a survival median of GBM patients who undergo the first-line standard treatment (surgery followed by adjuvant chemotherapy and local radiation) of 14 months since diagnosis, and a 5-year survival rate of only 6.8% [2, 3]. In addition, it is the most frequent of glioma tumors, accounting for 17% of this type of cancers [4]. Clinicians and researchers associate this high aggressiveness with its heterogeneity, rapid and dynamic progression and high invasive capacity [5, 6]. It is therefore clear that the characterization of GBM behavior is crucial for the development of therapeutic strategies against this cancer [7, 8].

The complexity of GBM evolution (and of cell biological processes in general), strongly depends on the particular microenvironment [9]. This makes it difficult the use of two-dimensional in vitro experiments (Petri dishes) to reproduce the actual behavior of cells in real tissues. In response to these limitations, microfluidics and micro-technologies permit to recreate the actual three-dimensional cell microenvironment in a much more realistic way, thus allowing to get results, associated with the complex biophysical and biochemical cell processes that govern the tumor dynamics, much closer to the actual in vivo situation than classical Petri dishes [10]. For instance, the study of tumor chemotaxis has been considerably developed [11, 12]. Consequently, and in general, these techniques allow testing drugs much more efficiently [13, 14].

Also, thanks to the flexibility, reproducibility, automation, integration and miniaturisation of microfluidic experiments, it is possible to generate great amounts of data, in the form of images and videos of cell culture evolution [15]. This opens the door to using Data Science and Machine Learning methods as new predicting tools [16, 17]. These tools can be considered therefore as a new paradigm in the analysis of complex, multifactorial, multiphysical and multiparametric biological problems, as those occurring during GBM evolution [18]. Therefore, it begins to be possible working in applied mathematical biology with complex nonlinear models, typical of other disciplines such as computational mechanics. In particular, models based on partial differential equations (PDEs) that have been used for predicting tumor evolution (see for instance [19]) may be enriched and adapted to consider more complex and coupled phenomena.

One example of particular interest is modeling of the so-called called go or grow paradigm [20], characteristic of many tumoral processes, among which GBM is a particular case. Indeed, hypoxia has been proposed as one of the main driving forces of GBM progression [21]. The fast proliferation of GBM cells close to brain blood vessels may provoke their occlusion, which in turns leads to local hypoxia. Thus, many cells die forming a necrotic core around the collapsed vessel. The surviving cells migrate towards other non collapsed vessels, looking for oxygenated areas. These migrating structures, formed of high-density areas of cells, are called pseudopalisades [22]. Once they reach a new vessel, their migration stops and the proliferation returns. This cyclic process of migration and proliferation is known in the scientific literature as the go or grow paradigm. It proposes that cells exhibit a migratory or proliferative activity depending on the oxygen concentration [20]. Hypoxia Inducible Factors (HIFs) are considered as the main biomolecular activators of this activity [23-25]. Recently, we have been able to reproduce these histological structures in vitro [26, 27]. Also, we developed a mathematical model incorporating the go or grow hypothesis, which allowed us to reproduce the GBM evolution under different experimental configurations also in vitro [28], and to derive, from cell culture images, information on the cell behavior [29].

However, parametric models are corseted by the mathematical relations that describe the a priori assumed biological hypotheses, so they present an obvious modeling bias. Besides, in our experiments, we try to understand the intrinsic mechanisms that control these biological processes; a knowledge that goes further than the numerical value of a specific parameter and is generally related to concepts with a clear biological meaning, such as whether the metabolic change is sharp or smooth, localized or distributed, or presents one or many different levels of transition.

A new promising family of neural networks, named as Physically-Guided Neural Networks with Internal Variables (PGNNIV), has just emerged as a tool to identify, evaluate and derive constitutive models from observable measurements [30, 31]. The fundamental idea is to incorporate the physical knowledge on the system into the network and to concentrate the learning power of Artificial Neural Networks in the intrinsic physical mechanisms that are intended to be found up. Very recently, a similar approach combining neural networks and physical equations (Physics-Informed Neural Networks (PINNs) [32]) has been proposed as a way to discover hidden mechanistic relationships using the Fisher–Kolmogorov–Petrovsky–Piskunov equation [33] as a benchmark problem, a concept that has been coined as Biologically-Informed Neural Networks (BINNs). All the same, PGNNIV offer greater flexibility in the definition of the internal variables of interest, including the non-measurable ones. Thus they are able to unravel more complex metabolic mechanisms, such as non-local or global models [31]. Additionally, they may deal with problems involving changeful external stimuli, that is, different source terms and boundary conditions, something that PINNs cannot afford.

In this work, we demonstrate how PGNNIV allow unravelling the mathematical structure that identifies the intrinsic metabolic mechanisms associated with cell changes due to the variation of some measurable fields, as the oxygen concentration around the cell that can be measured in microfluidic devices. This identification of the detailed go or grow mechanism related to hypoxia allows to accurately predict the cell evolution under highly variable external stimuli, including normoxic, gradient and hypoxic configurations, without the requirement of any initially prescribed parametric relation. The data used for testing the methodology are generated in silico, thus allowing evaluating the predictive and explanatory capacity of PGNNIV. Indeed, we demonstrate that this methodology does not only permit the identification of complex metabolic changes but also improves the prediction accuracy of parametric models.

The paper is structured as follows. The materials and methods section describes the formulation here presented. First, the mathematical model for Glioblastoma evolution is briefly revised, emphasising how the role of hypoxia is commonly taken into account for modeling the go or grow paradigm. Then, the model under the PGNNIV framework is presented, detailing the data generation and the network training process. In the Results section, the main results of the paper are presented: the unravelling of the metabolic behavior and the ability to predict the cell evolution under different oxygenation conditions. Finally, in the Discussion section, the present and future of this methodology is discussed, while in the Conclusions, the main results and conclusions are summarised.

Materials and methods

Our typical experimental configuration for the tests here analysed is shown in Fig 1. Here, the geometry of the model (assumed as one-dimensional, as the length of the lateral channels is much larger than the width of the chamber), and the different field variables are represented. Provided that the lengthscale is large enough, it is possible to identify the cell concentration with a continuum field u1 = u1(x, t). Besides, the oxygen concentration is associated with a field u2 = u2(x, t). The oxygen is supplied to the cell culture via the lateral channels. In response to this stimulus, cells will undergo migration and/or proliferation along the width of the chamber of length l.

Fig 1

Example of experimental configuration for modeling cell cultures.

Due to the long length of the lateral channels with respect to the width of the chamber, the geometry of the model is assumed as one-dimensional being the length the width of the chamber, l. The cell and the oxygen concentrations are associated with continuum fields u1 = u1(x, t) and u2 = u2(x, t) respectively. At the location of the lateral channels, x = 0, l, boundary conditions for both fields have to be supplied; we have represented a gradient configuration, but it is possible to feed oxygen by the two lateral channels in a symmetric configuration. (a) Scheme of the experimental configuration. (b) One-dimensional approximation of the cell culture. Created with BioRender.com.

Mathematical model of Glioblastoma cell culture evolution

Governing equations

Our starting point is a nonlinear reaction-diffusion system of partial differential equations that governs the evolution of GBM cells and the concentration of oxygen in a microfluidic device [28]. Although some works support the use of two [34-36] or even three [27] phenotypes for describing the cell population, we use a previously validated model offering some flexibility for modeling the switch between proliferative and migratory activity. In particular, the equations of the fields evolution are:

The first term of the R. H. S. of Eq (1a) represents the flow term associated with cell culture migration and has two contributions: the non-oriented motility term (modelled here as a random diffusion process) and the chemotaxis term , where cell motility is induced by the oxygen gradient. Πgo is a correction factor that will be discussed later. The second term corresponds to the reaction term and is associated with logistic growth [37], except for the correction term Πgr that will be also explained later.

With respect to the oxygen evolution equation, Eq (1b), the first term of the R. H. S. is again the flow term, consisting solely of oxygen diffusion, and the second corresponds to the oxygen consumption by cells. The correction between brackets in the second term is the Michaelis-Menten kinetic model [38] and accounts for the kinetics of oxidative phosphorylation that occurs in the membrane of cellular mitochondria [39].

Eqs (1a) and (1b) must be complemented with appropriate boundary and initial conditions. The initial condition is a known cell profile that is seeded in the microfluidic device at the beginning of the experiment (normally constant in the whole chamber). Note that we will refer to this time as t = 0 even if it is not necessarily the experiment starting time, identifying the instant when the cell culture profile is measured and the microfluidic device is fully oxygenated: with c(x) a given known function and the ambient oxygen level. The cell culture is subjected to a fixed oxygen concentration at the lateral channels. Besides, we assume that the walls of the culture chamber at the microfluidic devices are impermeable to cells, so only oxygen flow is allowed through them. In that case, the boundary conditions are: where we have defined as the cell flow, l is the length of the culture chamber and and are known functions defining the oxygen levels at the two lateral channels aside the chamber.

At this point, the presented framework has seven model parameters, D1, D2, χ, α1, α2, cs and km. Some of them have a well identified value in the scientific literature. For example:

The oxygen diffusion, D2 = 1 × 10−5 cm2 ⋅ s−1 has been reported in many works [40, 41].

The Michaelis-Menten constant, km = 2.5 mmHg, is very particular of the specific kinetics of the reaction in hands [41, 42].

All other parameters can be easily determined in specific well-controlled experiments. For example:

The parameters related to the logistic cell growth, α1 and cs, can be determined in cell growth experiments under fully oxygenated conditions and in absence of oxygen gradient, both in microfluidic devices [43, 44] or using cell spheroids [45].

The oxygen consumption rate, α2 is easily obtained by measuring the oxygen pressure at the ambient in an isolated system with a controlled cell culture population and for high oxygenation levels, such that the Michaelis-Menten correction, between brackets in Eq (1b), may be considered as 1. It is even possible to determine both km and α2 from the oxygen pressure using an Eadie–Hofstee diagram [46] or a Lineweaver–Burk plot [47].

The non-oxygen-mediated pedesis constant, D1, is more complicated to determine as spatial cell cultures are necessary. However, spheroid cultures [48] and microfluidic devices [26] offer a great opportunity for cell migration evaluation. If full oxygenation is guaranteed in the whole culture, no oxygen gradient is formed so non-oxygen mediated pedesis is easily computed, for instance, once given α1, D1 can be determined by evaluating the cell migration radial velocity V and using the Fisher’s model [49], .

The value of χ is substantially more difficult to determine. Indeed, as the cell migration depends on the oxygen level (and not only on the oxygen gradient), it is difficult to estimate this value from one single experiment or measurement. However, we can measure the cell culture migration under an oxygen gradient in a very localized region where the oxygen level may be considered almost constant [50]. Nevertheless, as it will be discussed later, we are rather interested in the overall value χΠgo. In this relation, χ is a reference value and Πgo is a correction term incorporating the effect of hypoxia in the migration.

Fig 2 illustrates some schematic experiments that can be implemented to determine the model parameters appearing in Eq (1) using conventional cell culture techniques and microfluidic devices.

Fig 2

Scheme of the different experiments that can be performed to obtain the model parameters.

Obtaining the source term parameters does not require a spatial cell distribution, although this is necessary for characterising the cell parameters associated with migration. (a) Determining α1 and cs. (b) Determining α2 and km. (c) Determining D1. (d) Determining χ. Created with BioRender.com.

In addition to the model parameters, the evolution of the GBM cell culture is also influenced by the boundary and initial conditions, in terms of the functions c(x), , x ∈ [0;l], and , , which play the role of problem data. That is, for the problem to be perfectly defined we need to specify the functions c, and together with the ambient oxygen pressure .

Go or grow activation functions

The metabolic behavior of the GBM cells, in particular, its response to changes in the oxygen pressure, is mathematically encoded in the functional form of Πgo and Πgr. These two functions regulate the activation/deactivation of both processes: migration and proliferation. There is sound evidence in the scientific literature that the switch between the proliferative and migratory activity in a cell population is hypoxia-mediated [23, 24], that is Πgo = Πgo(u2) and Πgr = Πgr(u2). However, there is not much knowledge about the details of this metabolic change (e.g.: are migration and proliferation simultaneous or not? Is the transition between them smooth? Is it strictly monotonic? Is it restricted to a narrow interval in the oxygen concentration region?).

In a recent paper [28], this transition was modelled by using piecewise linear functions of the ReLU kind, that is: and

Here, θgo and θgr play the role of oxygen thresholds. Additionally, it has been assumed that θgo = θgr, so this model implicitly assumes that Πgo(u2) + Πgr(u2) = 1, even if this consideration should, in principle, be modified to rely on a deeper understanding of the cell metabolism and in particular of the cell energy consumption. Indeed, the only biological evidence is that Πgr is a nondecreasing function and Πgo is a nonincreasing one. The relation between the two may be, in general, more complex, assumed as unknown, or even be unveiled also by the network.

The parameter values associated with Eqs (4) and (5) provided reasonably accurate results in the characterisation of certain cell cultures. In particular, GBM culture evolution of the cell line U251-MG in microfluidic devices has been well described, even for different experimental configurations using these expressions [28]. Similar results were obtained now using Machine Learning tools (in particular, using Convolutional Neural Networks), also revealing some limitations of the parametric model [29].

However, the go or grow model may differ from one GBM cell line to another. Besides, the model should be adapted for other tumor families or different frameworks. Therefore, since the functional relation u2 ↦ (Πgo, Πgr) encodes the cell metabolic changes in response to changes in the oxygen stimulus, its accurate characterisation is crucial for a complete understanding of the evolution of cell cultures, as it describes the changes that take place in the cell population behavior and consequently in the tumor progression [27, 51]. If Π = (Πgo, Πgr), the go or grow relation, may be written as: where is the unknown functional relation to be learned. Unravelling the one input-two output relation Π is therefore a key aspect in an in silico model able to capture tumor progression in an oxygenated medium. However, one main problem arises: as Πgo and Πgr are mathematical artefacts that only make sense when considered in Eq (1), they are non-measurable variables, so there is no experimental set up that permits to measure them directly. Furthermore, the measurement of the oxygen pressure in cell culture is usually difficult due to technical considerations, even if possible under some particular conditions [52, 53]. This adds an extra difficulty when defining or calibrating the model Π.

Physically-guided neural network with internal variables

Concept of PGNNIV

Physically-Guided Neural Networks with Internal Variables (PGNNIV) are a generalisation of the former concept of Physics Informed Neural Networks (PINN) [32]. In this latter, the physics of the problem informs the network via the output variables: the physical equations constrain the values of the output variables to belong to a certain physical manifold. For instance, to ensure that they satisfy a given partial differential equation. In other words, the loss function is directly defined in terms of the problem physics. PGNNIV go one step further, as in this case, the physical equations constrain the values reached by an arbitrary number of neurons in some intermediate layers. As a consequence, it is possible to interpret some hidden features and some relationships between internal variables (IV) of the problem that now acquire a physical interpretation [30, 31]. The physics does not constrain, but only guides the network learning capacity, as the measured data may be supplied to endow the network with explanatory capacity.

Going into the details, a PGNNIV is a problem formulated in the following archetypal way. Let us consider a PDE system of equations that is split into: where u and v are the unknown fields of the problem, and are functionals representing the known and unknown physical equations of the problem in hands. is a functional that specifies the boundary conditions, and f and g are known fields. Once discretised, Eq (7) has an analogous representation in finite-dimensional spaces in terms of vectorial functions , and and nodal values , , and . The Physically-Guided problem is therefore formulated as: where:

are the physical constraints, related to the relationships given by and .

and are functions specifying the input and the output of the problem, that is, the data used as starting point to make predictions and the data that we want to predict.

and are models. is the predictive model, whose aim is to infer accurate values for the output variables for a certain input set and is the explanatory model, whose objective is to unravel the hidden physics of the relation u ↦ v.

A PGNNIV is built when the problem (8) is formulated in the language of Neural Networks, with an appropriate structure and topology for and .

Discretised model

Spatial discretisation

Let us first discretise the Eq (1). This may be done by using Finite Differences (FD) or Finite Elements (FE) as it is usual when working with Partial Differential Equations (PDEs) [54]. The one-dimensional character and simple geometry of the cell culture in microfluidic devices under oxygen gradients [26] allow us to use FD to discretise the governing equations. Then, we define the nodal values of the fields u1 and u2 using the vectors 1 and 2, that is, u = u(x) where x = jΔx, j = 0, …, n, is the spatial coordinate associated with a given discretisation of the domain [0; l], . The spatial derivatives may be computed using any finite difference scheme, resulting in a linear operator . Then, Eq (1) results in:

We have used the symbols ⊙ and ⊘ for indicating pointwise multiplication and division respectively. It is important to note that Πgo and Πgr are here vector functions. The framework considered permits considering functional relationships, that is, the value of Π at a point x could depend on the value of the field u2 at the whole computational domain. However, the underlying nature of the go or grow framework allows us to consider Π(x) = Π(u2(x)), x ∈ [0; l], or equivalently, for the vector Π, Π = Π(u2), permitting to work with sparse graphs for the network topology (that is, sparse tensors and operators). In addition to the model Π, it would be possible to consider the rest of the specific model parameters (D1, χ, α1, cs, D2, α2, and km) as parameters to be learned during the training process. However, as this is similar to conventional parametric fitting (except for the fact that we use the broad capabilities of NN hardware and software [30]), we consider them here as known, with their values detailed next when specifying data generation.

In order to adapt the problem to our notations, let us describe Eq (9) as:

Temporal discretisation

With respect to the time integration, many options are possible. Multistep and Runge-Kutta methods [55] are one of the most efficient, although they are also computationally expensive. For our purposes, it is enough to consider a two-point scheme. Given the ODE , we discretise it using the generalized mid-point rule, that is, by approximating (t + Δt) − (t) ≃ Δt (β(t) + (1 − β)(t + Δt)), β ∈ [0; 1]. This approximation leads to the discretisation:

With this notation, taking β = 1 we recover the forward Euler method and with β = 0 we recover the backward Euler approach.

Applying Eq (11) to Eq (10) we obtain:

Finally, we may define the residual , that is indeed the equation encoding the problem physics, as:

The presented framework is generalisable to multistep and Runge-Kutta integrators. For instance, for the latter: with where a, b and c, i = 1, …, s are the particular coefficients of the selected numerical scheme.

In that case, the residual may be written as: where is given by Eq (15).

PGNNIV construction

The crucial part of building the PGNNIV is the definition of the network topology, as well as of the input and output layers. As explained when defining the biological problem, there is no way of straightforwardly measuring the variables Πgr and Πgo, so these will be our internal variables, = Π, while = (1, 2) are the measurable ones, that is, the cell and oxygen profiles at two consecutive time steps. Therefore:

The reason for defining the input and output variables this way is to achieve accurate predictive capacity, besides the required explanatory capacity. Indeed, once the model has been trained, it is possible to predict the outcome, that is, the cell and oxygen profiles at time t + Δt, from the ones given at time t. Note that the cell profiles (the output) at time t + Δt are obtained in real-time, as there is no need for solving any differential equation (we only need a network evaluation). The predictive and explanatory subnetworks are, therefore:

The PGNNIV graph and flow are illustrated in Fig 3. It is important to note that, although the explanatory subnetwork is the juxtaposition of two multilayer perceptrons, it is applied at each nodal value, so it acts as a convolutional network moving through the different collocation points in space. Note that if β = 1, as it is the case for this work, the input solely corresponds to the field values at time step n.

Fig 3

Structure of the PGNNIV.

Network topology and the different operators. The measurable (independent) variables, which are treated as the input variables, are represented in green while the predicted (dependent) variables, that are treated as the output of the network, are represented in magenta. The known operators are represented in blue, the unknown operators are represented in yellow and the hybrid operators are represented in orange. For illustrative purposes, we have split the network into three subnetworks. (a) Components of the network related to the problem physics, Eq (9), that is the discrete representation of Eq (1). (b) Explanatory subnetwork, whose aim is to unravel the relationship implicit in Eq (6). (c) Subnetwork related to the integration procedure, that is the one relating the input and output variables.

In a PGNNIV, the network loss function is the combination of a physics-associated term and a data-associated term. However, given the topology of the presented network, illustrated in Fig 3, all known physics of the problem is introduced explicitly in the network by means of the topology. Thus, the loss term is directly computed as: where we have denoted by the predicted value of the field by the PGNNIV. Recall that, observing the expression of the residual given by Eqs (13) and (19) may be written as: where we have defined the residual as a function of the data, = (, ).

Note that the architecture defining the PGNNIV, represented in Fig 3 and summarized in Eqs (18) and (20), representing respectively the predictive and explanatory network, does not depend on the initial and boundary conditions. This implies that the network may be trained with data coming from different experimental configurations, thus ensuring that the explanatory network is properly represented in the data. Once trained, predictions can be made, for any external stimuli and initial state. This is indeed one of the main differences between BINNs and PGNNIVs.

Data generation and training process

Data for model validation

Here we describe the data that will be used to feed the network. It is important to note that here the data-set is generated synthetically for validation purposes, but in real-life applications, this data-set would be the result of experimental measurements.

Benchmark models

In order to evaluate the performance of the method let us suppose four different true functional relationships for the metabolic model Π = Π(u2) that are described in Table 1. For illustration purposes, we assume for the different models that Πgr(x) = 1 − Πgo(x), although the PGNNIV may unravel the metabolic behavior for true models not satisfying this relationship.

Table 1

Different functional relationships defined for the validation procedure.

The different functions include features such as different smoothness and nonlinearities.

Chemotaxis activation function	Growth activation function	Parameters	Model name
Π_go(x)=I_[0;θ](x)	Π_gr(x)=I_[θ;+∞)(x)	θ	Heaviside, Binary step
Π_go(x)=(1-xθ)I_[0;θ](x)	Π_gr(x)=xθI_[0;θ](x)	θ	Piecewise linear, ReLU
Π_go(x)=kx+k	Π_gr(x)=xx+k	k	Michaelis-Menten
Π_go(x)=12(1+tanh(x-θΔθ))	Π_gr(x)=12(1-tanh(x-θΔθ))	θ, Δθ	Sigmoid, Logistic

We claim that our PGNNIV based on the governing Eq (1), encoding the known physics of the problem, is able to discover the actual biological metabolic model among the four presented in Table 1. This is possible due to the universal learning capabilities of neural networks [56, 57].

Profile generation

The data were generated by simulating cell profiles using Eq (1) with the boundary conditions (3). The model was first written using a dimensionless version, obtained by defining t = Tτ, x = Lξ, u1 = U1υ1 and u2 = U2υ2, where T is a characteristic time, L a characteristic length and U1 and U2 are characteristic cell and oxygen concentrations, obtaining: with boundary and initial conditions: where the dimensionless parameters and functions are:

In addition to the model parameters in Eq (24), we have to consider the ones related to the different go or grow models described in Table 1:

All parameters stated in Eqs (24) and (25) should have a precise biological meaning and depend on the problem physics. The values of the different parameters are reported in Table 2. Here, their value is only illustrative as they are used only for data generation, trying to make relevant all the biological phenomena.

Table 2

Dimensionless model parameters used for data generation.

The values are selected to make relevant all biological phenomena.

Parameter	Value
D¯_1	1
χ¯	1
α¯_1	0.5
c¯_s	10
D¯_2	1
α¯_2	0.05
k¯_m	2
l¯	50
θ¯	2
Δθ¯	2
k¯	2

The different cell and oxygen profiles were generated by using the method described in [58], especially suitable for parabolic partial differential equations. The system of equations was solved numerically by means of a time-space integrator based on a piecewise nonlinear Galerkin approach which is second order accurate in space, and compatible with this kind of nonlinear equations and boundary conditions, using the Matlab pdepe suit. Additional details may be found in [28]. A mesh size of Δξ = 1.0 and a time step of Δτ = 0.01 were used. As initial conditions, we set a value of and . The duration of the experiment is τ* = 10. Therefore, once the temporal series of the fields (cell and oxygen profile) are generated, the output of each simulation is an array of size [n, n, n] with n = τ*/Δτ + 1 = 1001, and n = 2.

Feeding the network

In order to recreate in silico different Glioblastoma On-Chip experiments, we created different cell profiles by varying the boundary conditions, that is, the oxygen levels and . Two families of configurations were simulated: symmetric and with oxygen gradient. The eleven in silico experiments performed are reported in Table 3. Each experiment is treated, from the PGNNIV point of view, as a batch of data as illustrated in Fig 3c: the batch k, k = 1, …, M, with M = 11, is therefore obtained by considering the k-th experimental configuration and the network is fed using each pair as input data and each pair as output, n = 0, …, n − 1. Each batch is therefore formed by an input of size [n − 1, n, n] (from values n = 0 to n = n − 2) and an output of size [n − 1, n, n] (from values n = 1 to n = n − 1). The objective is then learning the underlying go or grow model for a particular experimental condition.

Table 3

Experimental configurations used for data generation.

The different configurations recreate both symmetric and gradient configurations in low, medium and high oxygenated conditions (these values have to be compared with the model-associated ones, Eq (25)).

Configuration (k)	O¯2L	O¯2R
1	0	0
2	0	1
3	1	1
4	2	0
5	2	2
6	3	0
7	3	3
8	4	0
9	4	4
10	5	0
11	5	5

Training process

The neural network is trained using N = 103 epochs. At each epoch, all batches associated with the experimental configurations described in Table 3 are used for the network feeding: the M = 11 synthetic configurations, each one of them corresponding to one batch of data (so that the network parameters are updated after each batch feed), were sampled without replacement, so they were used in a different order at each epoch iteration. p = 80% of data at each batch is used for training purposes and 1 − p = 20% is used for testing the network. In total, N × M iterations of the network are considered until reaching convergence. The Adam optimizer [59] is selected with a learning rate r = 0.001 and an exponential decay rate of β1 = 0.8 for the first moment and β2 = 0.8 for the second moment are selected.

The training process takes 234 s in a core processor i7-6700 @3.4 GHz and RAM 64 GB, that is, without the use of GPU. Therefore, for an accurate comparison, we have adjusted the parametric model using also Adam optimizer, and it has taken 113 s.

Results

As in all problems involving PGNNIV, the neural network has both predictive and explanatory capacity. To illustrate both concepts, we will discuss first the explanatory capacity of the network, which is particular to this method. For comparative purposes, we will compare the learned relationship, in Eq (8) or Π in Fig 3 with standard parametric learning, where we postulate the model described by Eqs (4) and (5), also assuming θgo = θgr.

Then, we will comment on the predictive capacity of the network. As this capacity is not particular to the presented method but common to all regression techniques, we will compare our results to those obtained with standard parametric fitting.

Unravelling the metabolic changes of the GBM cells

In Fig 4 we depict the learned relationship Π for the four ground truth models proposed in Table 1. In all cases, a good agreement is shown between the real and the predicted models. Only when the parametric family for the go or grow model is adequately selected, the parametric learning (that is a particular PGNNIV where the function Π is parametrized) outperforms model-free PGNNIV, as it has been reported elsewhere [30, 31]. Note that when no explicit knowledge is assumed about cell metabolism, it is difficult to either derive or postulate parametric relations such as those in Eqs (4) and (5), which are solely used as a mere instrumental tool.

Fig 4

Unravelling capacity of the PGNNIV.

For all models in Table 1, presenting different nonlinearities, smoothness and scales, the ground truth model is recovered correctly, especially in relation to the growth metabolic behavior. (a) Heaviside transition model. (b) ReLU transition model. (c) Michaelis-Menten transition model. (d) Logistic transition model.

Fig 5 shows the errors when unravelling the metabolic behavior Π. Denoting by the model learned by the network, the error is defined as:

Fig 5

Error between the predicted and the real model.

The error of the prediction is robust over the different transition functions tested and outperforms any parametric fitting.

These errors are computed for both Πgo and Πgr. Except in the aforementioned case when the parametric family assumed includes the true model, the PGNNIV prediction clearly outperforms standard parametric approaches and keep the errors reduced for a broad family of families.

Predicting cell culture evolution

The aim now is to explore the predictive capacity of the neural network. Once the model Π has been learned, it is represented by the multilayer perceptron topology together with all the network parameters (weights and biases). Therefore it may be encapsulated as a one input—two output black box and inserted in any numerical integration scheme. For instance, we can consider any Runge-Kutta integrator of the form given by Eq (14) for the spatially-discretised equations, that is, to follow the approach for data generation, except for the fact that we use the learned model Π instead of any other of those presented in Table 1.

For illustrative purposes, let us compare the cell and oxygen profiles for three different boundary conditions: a normoxic configuration where , a hypoxic configuration where and a gradient configuration where and . Our aim is to predict the cell profile after τ = 20. The results are shown in Fig (6), where we have represented for each cell profile the real one (the derived when using directly the function in Table 1), the one predicted after fitting the parametric model described by Eqs (4) and (5) and the one predicted by PGNNIV. For all the models tested, a good agreement is shown between the predicted and the real profiles, and PGNNIV always outperforms the prediction of the parametric models, except, as it was explained before, for the ReLU case. The improvement of the prediction is particularly significant for the gradient configurations (Fig 6b). Indeed, the specific features of the model have a greater impact for oxygen levels in the transition between normoxic and hypoxic behavior since it is in this case when the differences between the different models most influence the cell evolution.

Fig 6

Prediction of the cell profile at τ* = 10 for the different models tested and different experimental configurations.

PGNNIV improves the prediction (when compared to the parametric model) of all cell profiles, slightly for the normoxic configuration and significantly for the gradient configuration. (a) Hypoxic configuration. (b) Gradient configuration. (c) Normoxic configuration.

In order to explore quantitatively the improvement, we define the error associated with the cell prediction as:

In Fig 7 we compare the error of the cell prediction given by Eq (27) when estimating the cell profile using the parametric approach and the one based on PGNNIV.

Fig 7

Prediction error.

The non-parametric model built using the PGNNIV approach is able to correctly estimate the cell profile evolution for any arbitrary model better than specific parametric ones. The predictions are more accurate (lower error) and more precise (lower variability) along with the different tested boundary conditions.

It is important to note that the training data-set was used for simulations for τ ≤ τ* with τ* = 10, so we explore here the prediction capacity of the network out of the region defined by the training data-set.

Computational requirements

In order to evaluate the computational requirements of the method, we compare it to standard parametric fitting performance. Recall that parametric fitting depends on the algorithm selected, in terms of both time and memory requirements. Conventional algorithms, such as Levenberg-Marquardt [60] are more memory demanding and require a large number of evaluations for large datasets, when compared for instance to Stochastic Gradient Descent (SGD) with small batch sizes. Therefore, for an accurate comparison, we have adjusted the parametric model using a similar PGNNIV (using Adam optimizer with same hyperparameters) except for the fact that the model function Π was totally parametrized in terms of θgo = θgr = θ.

The training process for the non-parametric PGNNIV takes 234 s in a core processor i7-6700 @3.4 GHz and RAM 64 GB, compared to the 113 s that takes for the parametric PGNNIV, using TensorFlow package (Python). The differences are obviously related to the number of network parameters, that are 2 × (8 × 8 + 8) = 144 for the non-parametric approach and 1 for the parametric one. Even so, in the two approaches, the low computational requirements are due to the NN framework used to face the problem.

Discussion and open possibilities

The present: Characterisation of complex biological cell processes

Discovery of hidden cell metabolisms is a major concern in biology. Indeed, unravelling the changes of the cell behavior when exposed to different stimuli can put us on the track of mechanisms driving the different cell signalling paths [61, 62]. As a result, parameters as the hypoxic threshold [28] are replaced by richer behaviors, as the ones illustrated in Fig 8.

Fig 8

Parametric vs non-parametric approaches.

The degree of information about the cell metabolism is richer in the non-parametric approach.

Moreover, from a mechanistic perspective, different metabolic paths and schemes may be tested in silico using computational approaches [63], in order to decide whether a path candidate is compatible with the metabolic changes discovered by means of the PGNNIV. Therefore, this knowledge on the macroscopic cell metabolic behavior at the population level is important not only from an epistemic point of view, for modeling purposes, but also as a promising tool for molecular biologists, in their attempt to isolate and define the different signalling pathways, thus providing a deeper understanding on the underlying mechanisms.

Sometimes, there are some energetic constraints (the more fundamental ones are those given by the first and second principles of thermodynamics) that restrict the accessible states in a biological system [64-66]. These constraints are translated into macroscopic ones in a continuum population model. For instance, one possible constraint is the former hypothesis that Πgr + Πgo = 1. However, this is a special case of the more general constraint G(Πgr, Πgo) = 0, that could be founded on an energetic argument about the resources available for the cell to grow or proliferate. All these extra constraints may be incorporated in the PGNNIV framework in a direct and straightforward fashion either by expressing some relational equations between variables explicitly, or by adding appropriate penalty terms in the loss function obliging to fulfil a mathematical constraint, such as p‖G(Πgr, Πgo)‖2, with p a penalty parameter [30].

A last remark is that PGNNIV, as any method inspired in Neural Networks, can be used as a universal approximator of the hidden interaction mechanisms between different cell populations, thus incorporating the ingredients of population sociology in systems biology [67]. For instance, if many cell populations are considered, C1, …, C, one may establish many ad-hoc interrelation mechanisms, λ = λ(C1, …, C), where λ is any model functional parameter, describing for instance migration or proliferation. The crosstalk between different cell populations has demonstrated to be important in many cellular processes such as those presented here [68, 69]. Of course, this interrelation would be properly learned if:

The known physics and biology of the system is well enough described in terms of specific mathematical equations. That is, all known mechanisms are explicitly stated, but only them. This enables the PGNNIV to concentrate its unravelling power in the unknown interrelations.

The available data is large enough to capture the prescribed dependency (as in the problem presented in this work). This is commonly difficult in many experimental sciences, particularly in biology, but new tools and trends such as microfluidics are promising in this regard.

As in any machine learning approach, care must be taken when interpreting the results and deriving conclusions, as the learning methods suffer from overfitting. A suitable strategy (train and test approaches, cross-validation, validation trials…) is therefore crucial for drawing generalisable conclusions.

The future: Towards in silico personalised medicine

This work presents a method for going from cell expression at the tissue level, that is, the formation of cellular structures such as pseudopalisades in GBM invasion, to cell behavior at the population level in response to the ambient stimuli. In one sense, it presents a link between clinicians and molecular biologists. A tumor biopsy may be extracted from one patient, cultivated and monitored in microfluidic devices where it may be exposed to different stimuli. The microfluidic devices have demonstrated to be capable of in silico reproducing histological features such as necrotic cores and pseudopalisades [26, 27], which can be captured using image and video techniques. PGNNIV is a tool able to infer, from the culture response to several stimuli, the intrinsic model of the cell reaction to such stimuli. This ability to integrate the knowledge of the response to different stimuli is a particular capability of PGNNIV, which puts them one step ahead of other methods such as BINNs in unravelling cell metabolism.

Once the intrinsic model is learned, its generalisable capability is much wider than the one offered by the histological features, as the latter is the response to very specific conditions (in mathematical terms, to very specific boundary conditions). In a sense, this strategy is the same as parametric fitting, except for the fact that there is no need for making any assumption about the model functional structure, provided that we know the internal variables we want to relate. This last issue is not minor, but is, indeed, the main objective of biology research: to make scientific conclusions about the effect and association of chemical factors with biological response.

This extra generalisation capability, when combined with appropriate mathematical models, offers new possibilities in in silico drug and treatment design. Once the cell response to the different stimuli is unravelled, we may variate the different stimuli levels and the different conditions. Again, in mathematical terms, this is represented by boundary conditions, initial conditions and source terms of the associated partial differential equations. Consequently, we can evaluate, again from a tissue perspective, more than from a cell population point of view, the general histological features that happen in the cell virtual tissue, which is called a Virtual Digital Twin. These histological features are responsible of key factors in cancer progression such as vessel occlusion [70], intravasation, extravasation and metastasis [70, 71], necrosis and activation of inflammatory response and/or angiogenesis [72, 73], in particular for GBM [74]. Therefore, in silico tests will accelerate the design of new drugs and therapies as they allow to test the hypotheses in a more flexible, faster and cheaper way.

It is important to note that all the here described sequence of steps relies on one single patient-specific biopsy, thus turning all this approach fully patient-specific, grounding this approach within the global framework of personalised in silico medicine [75]. From a specific patient, we infer specific histological and cellular features and, therefore, the different in silico trials are adapted to their particular disease, running away from parametric models, whose universality is corseted by their particular functional form.

The whole process is illustrated in Fig 9, where all the steps are summarized both in conceptual terms and using the particular example illustrated in this work. To conclude, PGNNIV is a tool for transferring the knowledge on the tissue response to knowledge on the cell gene expression. This knowledge is the one that enables to work on the cell expression level and to generalise both to arbitrary conditions and individual patients.

Fig 9

Summary of the described framework.

Starting from the clinical patient, we can recreate histological features associated with the tumoral tissue in the microfluidic devices, extract the cell response to stimuli using PGNNIV and use the information to evaluate different therapies and drug candidates. Created with BioRender.com.

Conclusions

The emerging technology of microfluidics has brought to the field of cell culture not only the possibilities of biotechnology research in more realistic biomimetic environments but also the chance of incorporating data-intensive tools such as Artificial Intelligence and Machine Learning. The term coined for this framework is Intelligent Microfluidics.

Here, we have illustrated how PGNNIV may be a valuable tool to infer the cell response at the population level from the cell response at the tumor microenvironment level in response to external stimuli. Unlike other Physically-Informed Data Science methods, in PGNNIV the physics does not constrain, but only guides the network learning capacity, as the measured data may be supplied to endow the network with explanatory capacity. Using in silico data, we have proven the unravelling capacity of PGNNIV for different benchmark test models, thus allowing us to work in model-free and non-parametric frameworks. Indeed, PGNNIV based simulations outperform the explanatory ability of any parametric model, except, at most, if the selected model belongs to the selected parametric family, which is, in practice, a strong and unrealistic assumption.

In addition, this explanatory ability is directly translated into an improvement of the prediction for different ambient conditions. The predictive ability does not only improve the one associated with classical parametric fitting but is also independent of the underlying model, thus making unnecessary the assumption of extra hypotheses, for example about cell metabolism. This improvement is achieved regardless of the environmental conditions, even if these ones are not used during the training process. In a certain sense, the ability of PGNNIV to predict out of the range of the training dataset is exploited here at its best.

The flexibility of PGNNIV allows that any information about biological systems may be incorporated totally or partially into the computations. This includes cell-cell interactions or cell-substrate interaction among other cues. In an opposite way, it is possible to focus the learning power on any relationship between fields of interest (biological, chemical, mechanical…) that is intended to be learned or quantified, both for theoretical (learn about the molecular metabolic processes) or for experimental (describe the main features of the process) purposes.

The presented methodology let us glimpse some steps in the direction of personalised medicine, as it is model-free, it allows to work with tissues extracted from different patients without the need of the specification of any particular model that would ruin out the method generalisation ability. Once characterised, the tumoral population may be in silico subjected to different stimuli and conditions, corresponding to different exploratory treatments. The response may be analysed from a clinician point of view, that is, at the tissue level, in order to evaluate the treatment success or failure. This strategy of in silico test-evaluation is quick and cheap and can strongly reduce animal experimentation, therefore facilitating research in areas such as biotechnology and biomedical engineering.

17 Dec 2021

Dear Dr Doblare,

Thank you very much for submitting your manuscript "Understanding glioblastoma invasion using the unravelling power of physically-guided neural networks with internal variables" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Inna Lavrik

Associate Editor

PLOS Computational Biology

Feilim Mac Gabhann

Editor-in-Chief

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The paper has a very well explain new method where it will help to modelling microfliudic techniques. The method is well explain with all the details and simulation results.

I found the paper well estructurare and written.

Reviewer #2: The review is uploaded as an attachment.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: No: The data and code is not available

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Pilar Guerrero

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

Submitted filename: Review.pdf

Click here for additional data file.

21 Feb 2022

Submitted filename: Letter_To_Reviewers.pdf

Click here for additional data file.

15 Mar 2022

Dear Dr Doblare,

We are pleased to inform you that your manuscript 'Understanding glioblastoma invasion using physically-guided neural networks with internal variables' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Inna Lavrik

Associate Editor

PLOS Computational Biology

Feilim Mac Gabhann

Editor-in-Chief

PLOS Computational Biology

***********************************************************

30 Mar 2022

PCOMPBIOL-D-21-01502R1

Understanding glioblastoma invasion using physically-guided neural networks with internal variables

Dear Dr Doblare,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Anita Estes

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol