Modeling and Simulation of VEGF Receptors Recruitment in Angiogenesis

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 themembrane stem from a rigorous thermodynamic setting, whereby strong simplifying assumptions are here taken and discussed. The multiphysics model is validated against experimental investigations.


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 membranebound reactants [1][2][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 diffusionreaction 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.

Mathematical Problems in Engineering
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.iiopen source finite element library, has been validated against codesigned experiments partially discussed in the companion paper [7].

Modeling VEGFR2 Diffusion Driven by
Its Specific Ligand

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 (R) and ligands (L) is described as a chemical reaction, which produces a receptorligand complex (C) where   and   are the kinetic constants of the forward and backward reaction, respectively.The reaction rate  (1) , measured in [/ 3 ], 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  → ℎ  ,  (1) , 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, div 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.

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][12][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), H[−] is the Heaviside step function,   = 72 ligands/m 2 is the concentration of substrate-immobilized ligand available for reaction (1),   is the time required for the complete mechanical deformation of the cell, V = ℓ/2  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  < /V; then, in the time span between  = /V and  = /V+, it increases rapidly from zero to   .We assume Figure 1: Surrogated mechanics: the cell-substrate contact dynamics is simulated by assuming that it is the substrate that gets deformed by the cell membrane, thus inducing a supply of ligands captured by function   in (5).
since complex is provided by  (1) only, and receptors are not generated.

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 as with  →  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 [0,   ] as where with  = {  ,   ,   } 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 V [0,  ] , V.However, the identification of these spaces falls beyond the scope of the present paper.

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

Energy Balance.
Denote with Ω the membrane, i.e., the spatial domain of problem.Consider an arbitrary region P ⊂ Ω.The first law of thermodynamics represents the balance of the interplay among the internal energy of P, the heat transferred in P and the power due to mass exchanged on P. The energy balance for the problem at hand reads dU d where Q  is the power due to heat transfer and T  is the power due to mass transfer.Denoting with P the bounding closed curve of P, and making reference to Figure 2 for the notation, they read The time variation of net internal energy U thus corresponds to the power expenditure of two external agents: a heat contribution Q  where   is the heat supplied by external agents and  →  is the heat flux vector; a mass contribution T  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) to The first law of thermodynamics is thus stated as follows: It must hold for any region P, since the latter is arbitrary.After simple algebra, the local form of the first principle thus reads where  is the net internal entropy of P,   is the entropy produced inside P,   is the entropy per unit time due to heat transfer, and   is the entropy per unit time due to mass transfer.The individual contributions read The 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 that Analogously 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 to By multiplying per  ≥ 0, replacing −  +div Ω [  →  ] by means of the energy balance (18), and some simple algebra, the entropy imbalance becomes Denote with the symbol   the quantity and with the symbol  (1) the following: By noting that one finally writes the entropy balance as where  →  =  →  + 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  (1) 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: (1)  (1) ≤ 0. (34b) 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) as A 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 u |  > 0 represents the average velocity of receptors when acted upon by a force of 1N/mol independent of the origin of the force.Definition (39) represents the physical requirement that both the pure (  = 0) and the saturated (  =    ) phases have vanishing mobilities.Neither the mobility u |  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 D |  = u |   is the receptor diffusivity.

Chemical Kinetics. The chemical kinetics of reaction (1)
is modeled via the law of mass action: At chemical equilibrium, as  (1) = 0 and  (1) = 0, the concentrations obey the relation which defines the constant of equilibrium  (1)  eq of reaction (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  (1) = 0, i.e., from ( 25) and ( 37) where Δ 0 =  0  −  0  −  0  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 becomes By subtracting (7b) from (7a), it comes out that with   from (5).These are the governing equations for   ,   , whereas   can be recovered from (45).Initial conditions read 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 as We make recourse to the Backward Euler method that leads to the following nonlinear problem in   (+Δ) and   (+Δ):

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.iiframework (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  = 1 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 ℓ = 20m was calculated from the measure of radius of 50 endothelial cells (ECs); receptor diffusivity D |  = 0.198m 2 /s 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  0  = 4.8molecules/m 2 .The kinetic parameters ligand/receptor were measured by surface plasmon resonance (SPR) measurements and the equilibrium constant turned out to be  (1)  eq = 354059.Finally, it has been taken    = 16000molecules/m 2 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.iiframework, we used it for the problem at hand.The mesh at time  = 85 is depicted in Figure 3(b).It consists of 14058 quadrilaterals, graded around the surface of contact between the cell and the substrate.
Figure 4 depicts the resulting evolution of the concentration of receptors in space and time.Figure 4(a) depicts a snapshot of   at  = 85 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  0  = 4.8receptors/m 2 .After 60 s, the concentration profile is perturbed at the bottom of the cell ( = 0), and such a perturbation "propagates" with time.This event can be clearly seen in Figure 4(a) at  = 85.At  = 600, 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 0.49receptors/m 2 and a steady state has not yet been reached.
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  = 7200 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

Normalized mass [∝ units]
Figure 5: Time evolution of the overall amount of bound VEGFR2gremlin complex on the membrane compared with fluorescence intensity from of the overall VEGFR2 (free and bound) in contact the substrate.The continuous lines correspond to the numerical simulation data, while points and error bars refer to the experimental data.To allow the comparison, both sets of data have been normalized to the respective values they reach at the final time  = 7200 s.The depletion of free receptors is clearly visible, whereas the sum of complex and free receptors concentration remains constant.
Mass is therefore conserved in our numerical validation.
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 ratecontrolling 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.

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.
) 2.4.2.Entropy Balance Equations.The second law of thermodynamics represents the balance of the interplay among the internal entropy of P and the entropy transferred in P due to mass exchange and heat transferred on P. The entropy balance for the problem at hand reads d d (P) − d  d (P) = Q  (P) + T  (P) ,

Figure 4 :
Figure 4: Time evolution of the spatial concentration   of free VEGFR2 along the cell membrane.(a) Snapshot at  = 85 s on a half of the membrane, with identification of points A, B, C, D, and E. (b) One-dimensional description of   at different times.Each curve plots the distribution of free receptors at different times  = 60, with  = 0, 1, 2, . ..,120 s from the beginning of the experiment at  = 0 to the final time  = 7200 s.Points A, B, C, D, and E correspond to (a).