Abstract
Angiogenesis, the process of new blood vessel formation from preexisting ones, plays a pivotal role in tumor growth. Vascular endothelial growth factor receptor-2 (VEGFR2) is the main proangiogenic tyrosine kinase receptor expressed by endothelial cells (ECs). VEGFR2 binds different ligands triggering vascular permeability and growth. VEGFR2-ligands accumulate in the extracellular matrix (ECM) and induce the polarization of ECs as well as the relocation of VEGFR2 in the basal cell membrane in contact with ECM. We propose here a multiphysical model to describe the dynamic of VEGFR2 on the plasma membrane. The governing equations for the relocation of VEGFR2 on the membrane stem from a rigorous thermodynamic setting, whereby strong simplifying assumptions are here taken and discussed. The multiphysics model is validated against experimental investigations.
1. Introduction
Vascular Endothelial Growth Factor Receptor-2 (VEGFR2) is a proangiogenic receptor expressed on endothelial cells (ECs) and is the main mediator of the angiogenic response. The interaction between VEGFR2 and extracellular ligands, produced by tumor cells, is essential to cancer growth. Specifically, ligand stimulation causes the relocation of VEGFR2 in the basal aspect in cells plated on ligand-enriched extracellular matrix both in vitro and in vivo, and ultimately receptors-ligands interaction activates the ECs division and proliferation towards tumor cells. Upon release, growth factors associate with the extracellular matrix and act as ECs guidance during neo-vessel formation.
Receptor-ligand interactions have been extensively studied and mathematical models have been proposed. Some concerned the estimation of the reaction rates for membrane-bound reactants [1–3]; a few models with different level of complexity account for adhesion receptor-ligand (as integrins and fibronectin) kinetics, receptor-ligand densities, cell rheology, and cytoskeletal force generation [4]. Only a few investigations concerned specifically VEGFR2 [5, 6].
Codesigned experiments and simulations for VEGFR2 have been recently developed, with biological findings and the predictive ability of the model extensively discussed in [7]. Here, we profoundly describe the modeling of VEGFR2 recruitment in angiogenesis, detailing the thermodynamic description, the weak formulation, and the algorithms for the numerical solution.
The mathematical model here proposed accounts for diffusion of VEGFR2 along the cellular membrane and for ligands-receptors chemical reactions. It is framed in the mechanics and thermodynamics of continua, following a general description proposed in [8], and takes advantage of successful descriptions of physically similar systems [9, 10]. The effect of the cell deformation on the diffusion-reaction process on the membrane is here strongly simplified, surrogating the effects of the change in geometry on the chemodiffusive equations with a fictitious source term of ligands, detailed in Section 2.2.
The model stems from continuity equations (for mass, energy, and entropy; see Section 2.1), standard chemical kinetics, summarized in Section 2.6, thermodynamic restrictions, and constitutive specifications, detailed in Section 2.4. This sequence provides the governing equations in a strong form in Section 2.7, which is converted in a weak form in Section 2.8 prior to the numerical approximation via the Finite Element Method (FEM). The partial differential equations of the model have been implemented in a computer code, with the ultimate goal to predict conditions for angiogenesis. The FEM code, implemented in the deal.ii open source finite element library, has been validated against codesigned experiments partially discussed in the companion paper [7].
2. Modeling VEGFR2 Diffusion Driven by Its Specific Ligand
2.1. Mass Balance Equations
A general formulation for the chemo-transport-mechanics problem is here tailored to model the relocation of VEGFR-2 driven by its specific ligand on the lipid bilayer membrane (henceforth denoted with ). The interaction between receptors and ligands is described as a chemical reaction, which produces a receptor-ligand complex where and are the kinetic constants of the forward and backward reaction, respectively. The reaction rate , measured in , quantifies the net formation of (C) as the difference between the forward and backward reaction rates.
Complex internalization and its return back to the surface are not considered in this model. Therefore, the mass balance equations are defined on the membrane as follows:Symbols in (2a)–(2c) have the following meaning (concentrations are defined in space and time, i.e., . The same holds for , , and . Functional dependence is specified when necessary only, to favor readability): (with ) is the molarity (i.e., the number of moles per unit area) of a generic species ; is the mass flux in terms of molecules, i.e., the number of molecules of species measured per unit length per unit time, and is a tangent vector field on the membrane; is the rate in number of molecules per unit volume per unit time at which species is generated by sources, and is the time.
Ligands, whose degradation is negligible, are immobilized in the substrate as they are in vitro. The complex are assumed to be immobile as well, i.e.,Since receptors are free to move along the membrane, reaction (1) portrays a conversion of mobile to trapped receptors and vice versa.
Equations (2a)–(2c) are defined on the cell membrane. Accordingly, the divergence operator has to be defined on the same surface. Denoting with the cell membrane unit normal,
Mass balance equations (2a)–(2c) shall be accompanied by the balance of force in order to model the mechanical deformation of the cell, whose boundary, the membrane, is the geometrical support of (2a)–(2c). Modeling the evolution of the Laplace-Beltrami operator that presides formulation concurrently with the large deformation of the cell is a phenomenally ambitious task, which is in progress motivated by the promising outcomes here shown. In the present work, we surrogate the mechanics with some simplifying assumptions.
2.2. Surrogated Mechanics
During the codesigned experimental test, the cell progressively spreads out on the substrate. Since the latter is enriched with immobilized ligands, the cell surface in contact with the support increases with time and results in a supply of available ligands for the chemical reaction (1) to occur. Mechanical models for cell spreading involve very sophisticated descriptions of active and passive behavior of cells [11–13], leading to simulations of impressive computational burden. In the present, seminal works do not account explicitly for the mechanical evolution of the cell, which keeps its original shape. Rather, we surrogate the effects of its change in geometry on the chemodiffusive equations (2a)–(2c) by introducing a source term of ligands whose expression is calibrated from experimental evidence [14]. The following expression for in (2b) is taken [7]:The path of reasoning beyond (5) is equivalent to consider the cell as rigid and the substrate much more deformable, so that the latter envelopes the spherical cell, as depicted in Figure 1.
In (5), is the Heaviside step function, ligands/μm2 is the concentration of substrate-immobilized ligand available for reaction (1), is the time required for the complete mechanical deformation of the cell, is the velocity of mechanical deformation (assumed to be constant until ), is the cell radius, is a parameter that identifies a finite time required for binding, is the curvilinear abscissa of our simplified geometry, and is the generic time. In view of (5), the supply of ligands at point on the membrane remains zero until ; then, in the time span between and , it increases rapidly from zero to . We assumesince complex is provided by only, and receptors are not generated.
In view of the above, mass balance equations (2a)–(2c) finally become
2.3. Weak Form and Boundary Conditions
The weak formulation of balance equations (7a)–(7c) comes out after multiplication by a suitable set of time independent test functions, here denoted with a superposed caret, and from an integration upon the membrane, exploiting Green’s formula to reduce the order of differentiation. Consider the mass balance (7a) as a prototype:In the former identity, a surface gradient operator arises in view of the integration by parts of the divergence term. Such a surface gradient, on the spherical smooth surface of the membrane, is defined aswith the cell membrane unit normal. Within weak formulations a contribution is usually defined on the boundary in view of the two-dimensional version of the divergence theorem. This is not the case for the cell membrane since it is a closed surface. The weak form of equations ((7b), (7c)) can be easily derived following the same path of reasoning.
In conclusion, the weak form of the balance equations can be written in the time interval aswherewith and . Column collects the time-dependent unknown fields. Column collects the steady-state test functions that correspond to the unknown fields in .
To computationally solve the (either weak (10) or strong (7a)–(7c)) problem, constitutive equations must be specified, which is the subject of Section 2.4. Ellipticity of the operators and functional and numerical properties of the solution and of its approximation depend on the constitutive assumptions and on the choice of the correct functional spaces . However, the identification of these spaces falls beyond the scope of the present paper.
2.4. Thermodynamics
In view of the assumptions made on the geometrical evolution of the membrane, there is no need to distinguish between material and spatial time derivative. When dealing with composite functions of the form , we will identify the total derivative with the roman symbol and the partial derivative with the symbol . It thus holds This notation will be used in the time derivative of internal and Helmholtz free energies, and of entropy.
2.4.1. Energy Balance
Denote with the membrane, i.e., the spatial domain of problem. Consider an arbitrary region . The first law of thermodynamics represents the balance of the interplay among the internal energy of , the heat transferred in and the power due to mass exchanged on . The energy balance for the problem at hand readswhere is the power due to heat transfer and is the power due to mass transfer. Denoting with the bounding closed curve of , and making reference to Figure 2 for the notation, they readThe time variation of net internal energy thus corresponds to the power expenditure of two external agents: a heat contribution where is the heat supplied by external agents and is the heat flux vector; a mass contribution in which the scalar denotes the change in specific energy provided by a unit supply of moles of species .
Since the geometry remains unchanged, one can define specific internal energy per unit mass or per unit surface, since none of them changes during the process. We choose to define it per unit surface, namely,Standard application of the surface divergence theorem and of mass balances (2a)–(2c) leads from (14a) and (14b) toThe first law of thermodynamics is thus stated as follows:It must hold for any region , since the latter is arbitrary. After simple algebra, the local form of the first principle thus reads
2.4.2. Entropy Balance Equations
The second law of thermodynamics represents the balance of the interplay among the internal entropy of and the entropy transferred in due to mass exchange and heat transferred on . The entropy balance for the problem at hand readswhere is the net internal entropy of , is the entropy produced inside , is the entropy per unit time due to heat transfer, and is the entropy per unit time due to mass transfer. The individual contributions readThe scalar denotes the change in specific entropy provided by a unit supply of moles of species . Equation (19) stems from the nontrivial assumption that mechanics does not contribute directly to the total entropy flow in the entropy balance equation. The second law of thermodynamics states thatAnalogously to the energy counterpart, we define the specific internal entropy per unit volume. Standard application of the divergence theorem and of mass balances (7a)–(7c) leads toBy multiplying per , replacing by means of the energy balance (18), and some simple algebra, the entropy imbalance becomesDenote with the symbol the quantityand with the symbol the following:By noting thatone finally writes the entropy balance as
2.4.3. Helmholtz Free Energy
The specific Helmholtz free energy is defined asand is taken as a function of temperature and concentrations, . It thus holdswhich can be inserted in (27) to derive the entropy imbalance in the final form:where
2.4.4. Thermodynamic Restrictions
Inequality (30) must hold for any region , since the latter was arbitrarily taken. Therefore, the following local inequality, usually termed after Clausius-Duhem, yields
This inequality must hold for any value of the time derivative of the temperature and of the concentrations , , and . Since they appear linearly in the inequality, the factors multiplying them must be zero, as otherwise it would be possible to find a value for the time derivatives that violate the inequality. Therefore, the following restrictions apply
In view of formula (32), the amount declared in (24) acquires the meaning of chemical potential and hence the term in (25) turns out to be the affinity of the reaction (1).
Equation (32) yields to the so-called Clausius-Plank inequality:that splits under the assumptions of Curie’s principle and thermal equilibrium in the following set of inequalities:
2.5. Constitutive Theory
We will assume henceforth that the system is in thermal equilibrium. The Helmholtz free energy density is furthermore additively decomposed into three separate parts:The free energy density of mobile guest atoms interacting with a host medium is described by an ideal solution model, which provides the following free energy density for the continuum approximation of mixing of the generic species :where is the ratio between the concentration and the saturation limit for each species. The chemical potential can be written accordingly to (32) asA strategy to satisfy the thermodynamic restriction (34a) is to model the flux of receptors by Fickian-diffusion, which linearly correlates to the gradient of its chemical potential :by means of a positive definite mobility tensor . The following isotropic nonlinear [15] specialization for the mobility tensor accounts for saturation. In formula (39), ; is the saturation limit for receptors. The mobility represents the average velocity of receptors when acted upon by a force of independent of the origin of the force. Definition (39) represents the physical requirement that both the pure () and the saturated () phases have vanishing mobilities. Neither the mobility nor the saturation concentration is assumed to change in time. Such a limitation can be removed without altering the conceptual picture if experimental data indicate an influence of temperature, stresses, or concentrations. Noting that Fick’s Law (38) specializes as follows:where is the receptor diffusivity.
2.6. Chemical Kinetics
The chemical kinetics of reaction (1) is modeled via the law of mass action:At chemical equilibrium, as and , the concentrations obey the relationwhich defines the constant of equilibrium of reaction (1).
2.6.1. Infinitely Fast Kinetics
Experimental evidences [7] show that (i) the equilibrium constant is high, thus favoring the formation of ligand-receptor complex and the depletion of receptors and ligands; (ii) the diffusion of receptors on the cell membrane is much slower than interaction kinetics. Accordingly, it can be assumed that the reaction kinetics is infinitely fast, in the sense that the time required to reach chemical equilibrium is orders of magnitude smaller than the time-scale of other processes. For these reasons, we assume that the concentrations of species are ruled by thermodynamic equilibrium at all times, and the concentration of complex is related to the others by the equation , i.e., from (25) and (37)where is the standard Gibbs free energy. Far from saturation, when ,having denoted with the following constant:
2.7. Governing Equations
Making use of (45) and (7a), (7c) finally becomesBy subtracting (7b) from (7a), it comes out thatwith from (5). These are the governing equations for , , whereas can be recovered from (45). Initial conditions read
2.8. Weak Form and Numerical Solution
The weak formulation in space results from multiplying the strong form of governing equations by a suitable set of tests functions and performing integration upon the domain. Specifically, the weak form (47), defining with a test function, readsApplying the divergence theorem over , the former equation transforms as follows:The weak form of (48), defining with a test function, reads after easy algebra:
The weak forms (51a) and (51b) can be transformed in a first-order Ordinary Differential Equation (ODE) in time if discretization is performed via separated variables, with spatial test and shape functions and nodal unknowns that depend solely on time. The usual Einstein summation convention is taken henceforth for repeated indexes.
The nonlinear ODEs read
For the time discretization of problems (53a) and (53b) finite difference schemes are generally used [16], for which the time derivatives of the concentrations are replaced by the finite differences asWe make recourse to the Backward Euler method that leads to the following nonlinear problem in and :
3. A Case Study
The chemo-transport model in its discrete form (55a) and (55b) was implemented in a finite element code with a fully coupled Newton-Raphson solver, first as a script in Wolfram Mathematica and later in the deal.ii framework ( http://www.dealii.org). The resulting code has been used to simulate the relocation of Vascular Endothelial Growth Factor Receptor-2 (a proangiogenic receptor expressed on endothelial cells, shortened in VEGFR2 henceforth) on the cell membrane during the mechanical adhesion of cells onto a ligand-enriched substrate, in a codesigned in vitro experiment detailed in [7].
Analyzing the evolution in time of VEGFR2 recruitment induced by immobilized ligands experimentally measured in [14] as fluorescence intensity of the overall VEGFR2 (free and bound), we calibrated = 600 s as the time for completion of the mechanical deformation and as the parameter that identifies a finite time required for binding. The remaining parameters for the simulation were calibrated by in vitro assays. The cell radius was calculated from the measure of radius of 50 endothelial cells (ECs); receptor diffusivity was obtained in [7] by Fluorescence Recovery after Photobleaching. The amount of VEGFR2 on cell membrane per area was calculated by dividing the number of high affinity binding sites [17, 18] for cell surface area and resulted in . The kinetic parameters ligand/receptor were measured by surface plasmon resonance (SPR) measurements and the equilibrium constant turned out to be . Finally, it has been taken and .
The initial discretization on the membrane is depicted in Figure 3(a). It consists of 6144 quadrilaterals, uniformly distributed on the spherical surface. The mesh is locally and automatically refined using the error indicator usually named after Kelly [19]. This error indicator tries to approximate the error per cell by integration of the jump of the gradient of the solution along the faces of each cell. Although developed for Laplace’s equation, Kelly’s error estimator has proven to be a suitable tool to generate locally refined meshes for a wide range of equations, not restricted to elliptic problems. For this reason and since a devoted class is available in the deal.ii framework, we used it for the problem at hand. The mesh at time is depicted in Figure 3(b). It consists of 14058 quadrilaterals, graded around the surface of contact between the cell and the substrate.
(a)
(b)
Figure 4 depicts the resulting evolution of the concentration of receptors in space and time. Figure 4(a) depicts a snapshot of at on the surface and locates points A-E on the membrane. Figure 4(b), which exploits the axis-symmetry of the problem, shows the concentration of free receptors on a diametric section. Concentrations range from zero up to the uniform initial concentration . After 60 s, the concentration profile is perturbed at the bottom of the cell (), and such a perturbation “propagates” with time. This event can be clearly seen in Figure 4(a) at . At , after the cell spreading, the basal portion of the cell membrane in contact with the substrate is essentially empty of free receptors because the equilibrium constant is very high and the reaction kinetics is assumed to be infinitely fast. Since no further supply is provided afterwards, the process becomes diffusion-dominated, and it slowly evolves towards a final steady state. At the end of the simulation, after two hours, the maximum concentration of free VEGFR2 is and a steady state has not yet been reached.
(a)
(b)
Immediately after the membrane gets in contact with the substrate, the evolution of is governed by chemistry. The increase of contact area, as results of cell deformation, further fuels reaction (1), until the cell spreading. Diffusion and reaction rule the evolution of receptor concentration afterwards. The time evolution of VEGFR2 recruitment in the basal portion of the cell membrane is depicted in Figure 5. It shows the overlay of the outcomes of simulations (continuous lines) and experiments (red dots displayed as intensity of fluorescence, which is proportional to the overall amount of bound and free receptors on the portion of membrane in contact with the substrate) of the dynamic recruitment of VEGFR2, normalized to the value it reaches at the final time s. Three phases of complex formation can be clearly identified in Figure 5 as time proceeds: (I) contact, (II) chemomechanical, and (III) diffusion. The depletion of free receptors is clearly visible, whereas the sum of complex and free receptors concentration must be a constant by definition. This numerical evidence is seen in Figure 5, showing that mass is conserved in our numerical validation.
The mechanical deformation of the cell is influenced by the chemical affinity of the VEGFR2-gremlin binding reaction coupled with intracellular actin dynamics. However, our current simplified model does not capture yet the complexity of such coupling. Instead, it surrogates the details of the mechanical deformation with the experimentally guided assumption of an empirical ad hoc time-sequence of gremlin supply on increasing portions of the cell membrane; see Section 2.2. This supply is faster than the contribution of diffusion, and the cell surface becomes depleted of receptors very rapidly where the cell adheres to the substrate. When the mechanical deformation terminates and the cell is eventually spread, the diffusion of receptors becomes the rate-controlling mechanism. During this final phase, receptors that diffuse across the boundary of the contact surface are immediately trapped and immobilized by the ligands on the substrate.
4. Conclusions
A continuum coupled model of transport-reaction has been dealt with in this paper. It describes the motion of receptors on the lipid membrane, with pointwise traps that account for the receptors-ligands-complex reaction. The model is framed in standard thermodynamics [20]. The energy and entropy contributions of the mass flux in the balance equations are accounted for. The selection of the Helmholtz free energy leads to the constitutive characterization, which coupled with balance equations provides the strong form of the problem of VEGF receptor recruitment in angiogenesis. The weak form of the problem is stated, eventually discretized in spatial finite elements with time discretization governed by the backward Euler method. Numerical simulations of the relocation of VEGFR2 on the cell membrane show the performances of the proposed approach, which in spite of its severe simplifying assumptions revealed predictive capability in a codesigned in vitro experiment detailed in [7].
This model allows identifying three different phases, biologically described for the first time, using cellular in vitro assays. Mathematical model clearly defines the temporal order of the chemical ligand-receptor interaction, the mechanical cell deformation, the receptor diffusion on cell membrane. Such mathematical method will be useful to predict the endothelial cell responses to different microenvironments, improving, e.g., the development of bioscaffold for tissue engineering.
Two strong assumptions are here introduced that will be removed in future works, namely, the hypothesis of infinitely fast kinetics (Section 2.6.1) and of rigid membrane with externally provided mass supply (Section 2.2). Indeed the positive outcomes achieved in capturing the membrane motion of receptors due to chemodiffusion with simpler surrogate models justify further efforts in adding model complexities, as complex internalization and cell spreading, to quote a few.
Conflicts of Interest
All the authors declare that no conflicts of interest exist.
Acknowledgments
V. Damioli acknowledges the support of Guido Berlucchi Foundation Young Researchers Mobility Programme. This work is supported by grants from Associazione Italiana per la Ricerca sul Cancro (IG 17276) to S. Mitola.