Skip to main content
Advertisement
  • Loading metrics

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

  • Jacobo Ayensa-Jiménez,

    Roles Conceptualization, Formal analysis, Investigation, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing

    Affiliations Mechanical Engineering Department, School of Engineering and Architecture, University of Zaragoza, Spain, Aragón Institute of Engineering Research (I3A), University of Zaragoza, Spain, Aragón Institute of Health Research (IIS Aragón), Spain

  • Mohamed H. Doweidar,

    Roles Funding acquisition, Supervision, Writing – review & editing

    Affiliations Mechanical Engineering Department, School of Engineering and Architecture, University of Zaragoza, Spain, Aragón Institute of Engineering Research (I3A), University of Zaragoza, Spain, Centro de Investigación Biomédica en Red en Bioingeniería, Biomateriales y Nanomedicina (CIBERBBN), Spain

  • Jose A. Sanz-Herrera,

    Roles Funding acquisition, Project administration, Resources, Supervision, Writing – review & editing

    Affiliation Mechanical Engineering Department, School of Engineering, University of Sevilla, Spain

  • Manuel Doblare

    Roles Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing – review & editing

    mdoblare@unizar.es

    Current address: Mariano Esquillor, S/N, Zaragoza, Spain

    Affiliations Mechanical Engineering Department, School of Engineering and Architecture, University of Zaragoza, Spain, Aragón Institute of Engineering Research (I3A), University of Zaragoza, Spain, Aragón Institute of Health Research (IIS Aragón), Spain, Centro de Investigación Biomédica en Red en Bioingeniería, Biomateriales y Nanomedicina (CIBERBBN), Spain

Abstract

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.

Author summary

In this work, we apply Physically-Guided Neural Networks with Internal Variables (PGNNIV) to the understanding of the Glioblastoma evolution process. We explain the metabolic changes between the proliferative and migrative activity of Glioblastoma cell cultures by using the go or grow activation functions as a pair of internal variables, whose dependence on the oxygen level is unravelled by some building blocks of the whole PGNNIV. Due to its model-free nature, our method is able to identify different classical mechanistic approaches and to outperform cell culture evolution predictions, as we demonstrate in the paper. Unlike Biologically-Informed Neural Networks we can assimilate data obtained from different boundary conditions and under different external stimuli to simulate the tumor progression under arbitrary conditions. We demonstrate this ability by comparing the predictions with different boundary conditions, resulting in different oxygenation conditions. This flexibility enables the use of our proposed method for personalised medical purposes, as the cell culture metabolic information, for a particular tumor, is encapsulated in a sub-network and may be used for arbitrary in silico tests.

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 [2325]. 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.

thumbnail
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.

https://doi.org/10.1371/journal.pcbi.1010019.g001

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 [3436] 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: (1a) (1b)

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: (2a) (2b) 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: (3a) (3b) (3c) (3d) 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.

thumbnail
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.

https://doi.org/10.1371/journal.pcbi.1010019.g002

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: (4) and (5)

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: (6) 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: (7) 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 F, G and H and nodal values u, v, f and g. The Physically-Guided problem is therefore formulated as: (8) where:

  • R are the physical constraints, related to the relationships given by F and G.
  • I and O 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 uv.

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 u1 and u2, that is, uij = ui(xj) where xj = 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 D. Then, Eq (1) results in: (9a) (9b)

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 Π, Πj = Πj(u2j), 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: (10)

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 y(t + Δt) − y(t) ≃ Δtf (βy(t) + (1 − β)y(t + Δt)), β ∈ [0; 1]. This approximation leads to the discretisation: (11)

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: (12)

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

The presented framework is generalisable to multistep and Runge-Kutta integrators. For instance, for the latter: (14) with (15) where aij, bi and ci, i = 1, …, s are the particular coefficients of the selected numerical scheme.

In that case, the residual R may be written as: (16) where ki 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, v = Π, while u = (u1, u2) are the measurable ones, that is, the cell and oxygen profiles at two consecutive time steps. Therefore: (17a) (17b)

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: (18a) (18b)

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.

thumbnail
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.

https://doi.org/10.1371/journal.pcbi.1010019.g003

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: (19) where we have denoted by the predicted value of the field u by the PGNNIV. Recall that, observing the expression of the residual given by Eqs (13) and (19) may be written as: (20) where we have defined the residual as a function of the data, R = R(un, un+1).

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.

thumbnail
Table 1. Different functional relationships defined for the validation procedure.

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

https://doi.org/10.1371/journal.pcbi.1010019.t001

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 = , x = , 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: (21a) (21b) with boundary and initial conditions: (22a) (22b) (23a) (23b) (23c) (23d) where the dimensionless parameters and functions are: (24)

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: (25)

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.

thumbnail
Table 2. Dimensionless model parameters used for data generation.

The values are selected to make relevant all biological phenomena.

https://doi.org/10.1371/journal.pcbi.1010019.t002

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 [nt, nx, nu] with nt = τ*/Δτ + 1 = 1001, and nu = 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, …, nt − 1. Each batch is therefore formed by an input of size [nt − 1, nx, nu] (from values n = 0 to n = nt − 2) and an output of size [nt − 1, nx, nu] (from values n = 1 to n = nt − 1). The objective is then learning the underlying go or grow model for a particular experimental condition.

thumbnail
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)).

https://doi.org/10.1371/journal.pcbi.1010019.t003

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, H 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.

thumbnail
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.

https://doi.org/10.1371/journal.pcbi.1010019.g004

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

thumbnail
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.

https://doi.org/10.1371/journal.pcbi.1010019.g005

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.

thumbnail
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.

https://doi.org/10.1371/journal.pcbi.1010019.g006

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

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.

thumbnail
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.

https://doi.org/10.1371/journal.pcbi.1010019.g007

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.

thumbnail
Fig 8. Parametric vs non-parametric approaches.

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

https://doi.org/10.1371/journal.pcbi.1010019.g008

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 [6466]. 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 Ggr, Π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 pGgr, Π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, …, Cn, one may establish many ad-hoc interrelation mechanisms, λ = λ(C1, …, Cn), 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.

thumbnail
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.

https://doi.org/10.1371/journal.pcbi.1010019.g009

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.

References

  1. 1. Organization WH, et al. WHO report on cancer: setting priorities, investing wisely and providing care for all. 2020;.
  2. 2. Ostrom QT, Gittleman H, Farah P, Ondracek A, Chen Y, Wolinsky Y, et al. CBTRUS statistical report: Primary brain and central nervous system tumors diagnosed in the United States in 2006-2010. Neuro-oncology. 2013;15(suppl_2):ii1–ii56. pmid:24137015
  3. 3. Oike T, Suzuki Y, Sugawara Ki, Shirai K, Noda Se, Tamaki T, et al. Radiotherapy plus concomitant adjuvant temozolomide for glioblastoma: Japanese mono-institutional results. PLoS One. 2013;8(11):e78943. pmid:24265731
  4. 4. Brat DJ. Glioblastoma: biology, genetics, and behavior. American Society of Clinical Oncology Educational Book. 2012;32(1):102–107. pmid:24451717
  5. 5. Friedmann-Morvinski D. Glioblastoma heterogeneity and cancer cell plasticity. Critical Reviews™ in Oncogenesis. 2014;19(5). pmid:25404148
  6. 6. Furnari FB, Fenton T, Bachoo RM, Mukasa A, Stommel JM, Stegh A, et al. Malignant astrocytic glioma: genetics, biology, and paths to treatment. Genes & development. 2007;21(21):2683–2710. pmid:17974913
  7. 7. Xie Q, Mittal S, Berens ME. Targeting adaptive glioblastoma: an overview of proliferation and invasion. Neuro-oncology. 2014;16(12):1575–1584. pmid:25082799
  8. 8. Kathagen-Buhmann A, Schulte A, Weller J, Holz M, Herold-Mende C, Glass R, et al. Glycolysis and the pentose phosphate pathway are differentially associated with the dichotomous regulation of glioblastoma cell migration versus proliferation. Neuro-oncology. 2016;18(9):1219–1229. pmid:26917237
  9. 9. Bray D. Cell movements: from molecules to motility. Garland Science; 2000.
  10. 10. Sackmann EK, Fulton AL, Beebe DJ. The present and future role of microfluidics in biomedical research. Nature. 2014;507(7491):181–189. pmid:24622198
  11. 11. Mosadegh B, Saadi W, Wang SJ, Jeon NL. Epidermal growth factor promotes breast cancer cell chemotaxis in CXCL12 gradients. Biotechnology and bioengineering. 2008;100(6):1205–1213. pmid:18553401
  12. 12. Tatárová Z, Abbuehl J, Maerkl S, Huelsken J. Microfluidic co-culture platform to quantify chemotaxis of primary stem cells. Lab on a Chip. 2016;16(10):1934–1945. pmid:27137768
  13. 13. Shin Y, Han S, Jeon JS, Yamamoto K, Zervantonakis IK, Sudo R, et al. Microfluidic assay for simultaneous culture of multiple cell types on surfaces or within hydrogels. Nature protocols. 2012;7(7):1247–1259. pmid:22678430
  14. 14. Guckenberger DJ, de Groot TE, Wan AM, Beebe DJ, Young EW. Micromilling: a method for ultra-rapid prototyping of plastic microfluidic devices. Lab on a Chip. 2015;15(11):2364–2378. pmid:25906246
  15. 15. Coluccio ML, Perozziello G, Malara N, Parrotta E, Zhang P, Gentile F, et al. Microfluidic platforms for cell cultures and investigations. Microelectronic Engineering. 2019;208:14–28.
  16. 16. Riordon J, Sovilj D, Sanner S, Sinton D, Young EW. Deep learning with microfluidics for biotechnology. Trends in biotechnology. 2019;37(3):310–324. pmid:30301571
  17. 17. Galan EA, Zhao H, Wang X, Dai Q, Huck WT, Ma S. Intelligent Microfluidics: The Convergence of Machine Learning and Microfluidics in Materials Science and Biomedicine. Matter. 2020;3(6):1893–1922.
  18. 18. Cai X, Briggs RG, Homburg HB, Young IM, Davis EJ, Lin YH, et al. Application of microfluidic devices for glioblastoma study: current status and future directions. Biomedical Microdevices. 2020;22(3):1–10. pmid:32870410
  19. 19. Eikenberry SE, Sankar T, Preul M, Kostelich E, Thalhauser C, Kuang Y. Virtual glioblastoma: growth, migration and treatment in a three-dimensional mathematical model. Cell proliferation. 2009;42(4):511–528. pmid:19489983
  20. 20. Hatzikirou H, Basanta D, Simon M, Schaller K, Deutsch A. ‘Go or grow’: the key to the emergence of invasion in tumour progression? Mathematical medicine and biology: a journal of the IMA. 2012;29(1):49–65. pmid:20610469
  21. 21. Brat DJ, Van Meir EG. Vaso-occlusive and prothrombotic mechanisms associated with tumor hypoxia, necrosis, and accelerated growth in glioblastoma. Laboratory Investigation. 2004;84(4):397. pmid:14990981
  22. 22. Brat DJ, Castellano-Sanchez AA, Hunter SB, Pecot M, Cohen C, Hammond EH, et al. Pseudopalisades in glioblastoma are hypoxic, express extracellular matrix proteases, and are formed by an actively migrating cell population. Cancer research. 2004;64(3):920–927. pmid:14871821
  23. 23. Lu X, Kang Y. Hypoxia and hypoxia-inducible factors: master regulators of metastasis. Clinical cancer research. 2010;16(24):5928–5935. pmid:20962028
  24. 24. Carreau A, Hafny-Rahbi BE, Matejuk A, Grillon C, Kieda C. Why is the partial oxygen pressure of human tissues a crucial parameter? Small molecules and hypoxia. Journal of cellular and molecular medicine. 2011;15(6):1239–1253. pmid:21251211
  25. 25. Wang P, Yan Q, Liao B, Zhao L, Xiong S, Wang J, et al. The HIF1α/HIF2α-miR210-3p network regulates glioblastoma cell proliferation, dedifferentiation and chemoresistance through EGF under hypoxic conditions. Cell death & disease. 2020;11(11):1–13. pmid:33208727
  26. 26. Ayuso JM, Virumbrales-Muñoz M, Lacueva A, Lanuza PM, Checa-Chavarria E, Botella P, et al. Development and characterization of a microfluidic model of the tumour microenvironment. Scientific reports. 2016;6(1):1–16. pmid:27796335
  27. 27. Ayuso JM, Monge R, Martínez-González A, Virumbrales-Muñoz M, Llamazares GA, Berganzo J, et al. Glioblastoma on a microfluidic chip: Generating pseudopalisades and enhancing aggressiveness through blood vessel obstruction events. Neuro-oncology. 2017;19(4):503–513. pmid:28062831
  28. 28. Ayensa-Jiménez J, Pérez-Aliacar M, Randelovic T, Oliván S, Fernández L, Sanz-Herrera JA, et al. Mathematical formulation and parametric analysis of in vitro cell models in microfluidic devices: application to different stages of glioblastoma evolution. Scientific Reports. 2020;10(1):1–21. pmid:33273574
  29. 29. Pérez-Aliacar M, Doweidar MH, Doblaré M, Ayensa-Jiménez J. Predicting cell behaviour parameters from glioblastoma on a chip images. A deep learning approach. Computers in Biology and Medicine. 2021; p. 104547. pmid:34139437
  30. 30. Ayensa-Jiménez J, Doweidar MH, Sanz-Herrera JA, Doblaré M. Identification of state functions by physically-guided neural networks with physically-meaningful internal layers. arXiv preprint arXiv:201108567. 2020;.
  31. 31. Ayensa-Jiménez J, Doweidar MH, Sanz-Herrera JA, Doblaré M. On the application of Physically-Guided Neural Networks with Internal Variables to Continuum Problems. arXiv preprint arXiv:201111376. 2020;.
  32. 32. Raissi M, Perdikaris P, Karniadakis GE. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics. 2019;378:686–707.
  33. 33. Lagergren JH, Nardini JT, Baker RE, Simpson MJ, Flores KB. Biologically-informed neural networks guide mechanistic modeling from sparse experimental data. PLoS computational biology. 2020;16(12):e1008462. pmid:33259472
  34. 34. Fedotov S, Iomin A. Migration and proliferation dichotomy in tumor-cell invasion. Physical Review Letters. 2007;98(11):118101. pmid:17501094
  35. 35. Gerlee P, Nelander S. The impact of phenotypic switching on glioblastoma growth and invasion. PLoS computational biology. 2012;8(6):e1002556. pmid:22719241
  36. 36. Stepien TL, Rutter EM, Kuang Y. Traveling waves of a go-or-grow model of glioma growth. SIAM Journal on Applied Mathematics. 2018;78(3):1778–1801.
  37. 37. Tsoularis A, Wallace J. Analysis of logistic growth models. Mathematical biosciences. 2002;179(1):21–55. pmid:12047920
  38. 38. Cornish-Bowden A. The origins of enzyme kinetics. FEBS letters. 2013;587(17):2725–2730. pmid:23791665
  39. 39. Tang PS. On the rate of oxygen consumption by tissues and lower organisms as a function of oxygen tension. The Quarterly Review of Biology. 1933;8(3):260–274.
  40. 40. Tannock IF. Oxygen diffusion and the distribution of cellular radiosensitivity in tumours. The British journal of radiology. 1972;45(535):515–524. pmid:5067983
  41. 41. Daşu A, Toma-Daşu I, Karlsson M. Theoretical simulation of tumour oxygenation and results from acute and chronic hypoxia. Physics in Medicine & Biology. 2003;48(17):2829. pmid:14516104
  42. 42. Alper T, Howard-Flanders P. Role of oxygen in modifying the radiosensitivity of E. coli B. Nature. 1956;178(4540):978–979. pmid:13378491
  43. 43. Lei KF, Wu MH, Hsu CW, Chen YD. Real-time and non-invasive impedimetric monitoring of cell proliferation and chemosensitivity in a perfusion 3D cell culture microfluidic chip. Biosensors and Bioelectronics. 2014;51:16–21. pmid:23920091
  44. 44. Tao FF, Xiao X, Lei KF, Lee IC. Based cell culture microfluidic system. BioChip Journal. 2015;9(2):97–104.
  45. 45. Vinci M, Gowan S, Boxall F, Patterson L, Zimmermann M, Lomas C, et al. Advances in establishment and analysis of three-dimensional tumor spheroid-based functional assays for target validation and drug evaluation. BMC biology. 2012;10(1):1–21. pmid:22439642
  46. 46. Hofstee B. Non-inverted versus inverted plots in enzyme kinetics. Nature. 1959;184(4695):1296–1298. pmid:14402470
  47. 47. Lineweaver H, Burk D. The determination of enzyme dissociation constants. Journal of the American chemical society. 1934;56(3):658–666.
  48. 48. Ayuso JM, Basheer HA, Monge R, Sánchez-Álvarez P, Doblaré M, Shnyder SD, et al. Study of the chemotactic response of multicellular spheroids in a microfluidic device. PloS one. 2015;10(10):e0139515. pmid:26444904
  49. 49. Fisher RA. The wave of advance of advantageous genes. Annals of eugenics. 1937;7(4):355–369.
  50. 50. Funamoto K, Zervantonakis IK, Liu Y, Ochs CJ, Kim C, Kamm RD. A novel microfluidic platform for high-resolution imaging of a three-dimensional cell culture under a controlled hypoxic environment. Lab on a chip. 2012;12(22):4855–4863. pmid:23023115
  51. 51. Monteiro AR, Hill R, Pilkington GJ, Madureira PA. The role of hypoxia in glioblastoma invasion. Cells. 2017;6(4):45. pmid:29165393
  52. 52. Lam SF, Shirure VS, Chu YE, Soetikno AG, George SC. Microfluidic device to attain high spatial and temporal control of oxygen. PLoS One. 2018;13(12):e0209574. pmid:30571786
  53. 53. Zirath H, Rothbauer M, Spitz S, Bachmann B, Jordan C, Müller B, et al. Every breath you take: non-invasive real-time oxygen biosensing in two-and three-dimensional microfluidic cell models. Frontiers in physiology. 2018;9:815. pmid:30018569
  54. 54. Ganesan S, Lingeshwaran S. Galerkin finite element method for cancer invasion mathematical model. Computers & Mathematics with Applications. 2017;73(12):2603–2617.
  55. 55. Lambert JD. Numerical methods for ordinary differential systems: the initial value problem. John Wiley & Sons, Inc.; 1991.
  56. 56. Cybenko G. Approximations by superpositions of a sigmoidal function. Mathematics of Control, Signals and Systems. 1989;2:183–192.
  57. 57. Hornik K. Approximation capabilities of multilayer feedforward networks. Neural networks. 1991;4(2):251–257.
  58. 58. Skeel RD, Berzins M. A method for the spatial discretization of parabolic equations in one space variable. SIAM journal on scientific and statistical computing. 1990;11(1):1–32.
  59. 59. Kingma DP, Ba J. Adam: A method for stochastic optimization. arXiv preprint arXiv:14126980. 2014;.
  60. 60. Levenberg K. A method for the solution of certain non-linear problems in least squares, Quarterly of applied mathematics. 1944; 2(2):164–168.
  61. 61. Lavrik I, Golks A, Krammer PH. Death receptor signaling. Journal of cell science. 2005;118(2):265–267. pmid:15654015
  62. 62. Laplante M, Sabatini DM. mTOR signaling at a glance. Journal of cell science. 2009;122(20):3589–3594. pmid:19812304
  63. 63. Kholodenko BN. Cell-signalling dynamics in time and space. Nature reviews Molecular cell biology. 2006;7(3):165–176. pmid:16482094
  64. 64. Benzinger T. Thermodynamics, chemical reactions and molecular biology. Nature. 1971;229(5280):100–102. pmid:4923089
  65. 65. Haynie DT. Biological thermodynamics. Cambridge University Press; 2001.
  66. 66. Dill K, Bromberg S. Molecular driving forces: statistical thermodynamics in biology, chemistry, physics, and nanoscience. Garland Science; 2010.
  67. 67. Ganesh S, Utebay B, Heit J, Coskun AF. Cellular sociology regulates the hierarchical spatial patterning and organization of cells in organisms. Open Biology. 2020;10(12):200300. pmid:33321061
  68. 68. Oliveira AI, Anjo SI, de Castro JV, Serra SC, Salgado AJ, Manadas B, et al. Crosstalk between glial and glioblastoma cells triggers the “go-or-grow” phenotype of tumor cells. Cell Communication and Signaling. 2017;15(1):1–12.
  69. 69. Chen JWE, Lumibao J, Leary S, Sarkaria JN, Steelman AJ, Gaskins HR, et al. Crosstalk between microglia and patient-derived glioblastoma cells inhibit invasion in a three-dimensional gelatin hydrogel model. Journal of neuroinflammation. 2020;17(1):1–15.
  70. 70. Weis S, Cui J, Barnes L, Cheresh D. Endothelial barrier disruption by VEGF-mediated Src activity potentiates tumor cell extravasation and metastasis. The Journal of cell biology. 2004;167(2):223–229. pmid:15504909
  71. 71. Chiang SP, Cabrera RM, Segall JE. Tumor cell intravasation. American Journal of Physiology-Cell Physiology. 2016;311(1):C1–C14. pmid:27076614
  72. 72. Carmeliet P, Jain RK. Angiogenesis in cancer and other diseases. nature. 2000;407(6801):249–257. pmid:11001068
  73. 73. Greten FR, Grivennikov SI. Inflammation and cancer: triggers, mechanisms, and consequences. Immunity. 2019;51(1):27–41. pmid:31315034
  74. 74. Rong Y, Durden DL, Van Meir EG, Brat DJ. ‘Pseudopalisading’necrosis in glioblastoma: a familiar morphologic feature that links vascular pathology, hypoxia, and angiogenesis. Journal of Neuropathology & Experimental Neurology. 2006;65(6):529–539.
  75. 75. Cirillo D, Valencia A. Big data analytics for personalized medicine. Current opinion in biotechnology. 2019;58:161–167. pmid:30965188