# Modeling and Simulation of VEGF Receptors Recruitment in Angiogenesis

Authors: A. Salvadori; V. Damioli; C. Ravelli; S. Mitola
Journal: Mathematical Problems in Engineering
 (2018)

---

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


---

## Body

## 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 Section2.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 (R) and ligands (L) is described as a chemical reaction, which produces a receptor-ligand complex (C)(1)R+L⇄kbkfC,where kf and kb are the kinetic constants of the forward and backward reaction, respectively. The reaction rate w(1), measured in [mol/m3s], 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:(2a)∂cR∂t+divΩ⁡h→R+w1=sR,(2b)∂cL∂t+divΩ⁡h→L+w1=sL,(2c)∂cC∂t+divΩ⁡h→C-w1=sC.Symbols in (2a)–(2c) have the following meaning (concentrations cβ are defined in space and time, i.e., cβ=cβx→,t. The same holds for h→β, w1, and sβ. Functional dependence is specified when necessary only, to favor readability): cβ (with β=R,L,C ) is the molarity (i.e., the number of moles per unit area) of a generic species β; h→β 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; sβ is the rate in number of molecules per unit volume per unit time at which species β is generated by sources, and t 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.,(3)h→L=h→C=0→.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 n→ the cell membrane unit normal,(4)divΩ⁡h→R=n→·curln→×h→R=div⁡h→R-∇h→Rn→·n→.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 sL whose expression is calibrated from experimental evidence [14]. The following expression for sL in (2b) is taken [7]:(5)sLx,t=cL¯t¯Ht-xvHt¯-t+xvThe 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.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 functionsL in (5).In (5), H[-] is the Heaviside step function, cL¯=72 ligands/μm2 is the concentration of substrate-immobilized ligand available for reaction (1), tf is the time required for the complete mechanical deformation of the cell, v=πl/2tf is the velocity of mechanical deformation (assumed to be constant until tf), l is the cell radius, t¯≪tf is a parameter that identifies a finite time required for binding, x is the curvilinear abscissa of our simplified geometry, and t is the generic time. In view of (5), the supply of ligands at point x on the membrane remains zero until t<x/v; then, in the time span between t=x/v and t=x/v+t¯, it increases rapidly from zero to cL¯. We assume(6)sR=sC=0,since complex is provided by w(1) only, and receptors are not generated.In view of the above, mass balance equations (2a)–(2c) finally become(7a)∂cR∂t+divΩ⁡h→R+w1=0(7b)∂cL∂t+w1-sLx,t=0(7c)∂cC∂t-w1=0


### 2.3. Weak Form and Boundary Conditions

The weak formulation of balance equations (7a)–(7c) comes out after multiplication by a suitable set oftime 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:(8)∫Ωc^R∂cR∂t+divΩ⁡h→R+w1dA=∫Ωc^R∂cR∂tdA+∫ΩdivΩc^Rh→R-∇Ωc^R·h→RdA+∫Ωc^Rw1dA=∫Ωc^R∂cR∂tdA-∫Ω∇Ωc^R·h→RdA+∫Ωc^Rw1dA=0.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(9)∇Ωc^R=∇c^R-n→·∇c^Rn→with n→ 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 interval0,tf as(10)Findy∈V0,tfsuch that∂∂tby^,yt+ay^,yt=fy^∀y^∈Vwhere(11)by^,y=∫Ωc^RcR+c^LcL+c^CcCdA,ay^,y=-∫Ω∇Ωc^R·h→RdA+∫Ωc^R+c^L-c^Cw1dA,fy^=∫Ωc^LsLdAwith y=cR,cL,cC and y^=c^R,c^L,c^C. Column y collects the time-dependent unknown fields. Column y^ collects the steady-state test functions that correspond to the unknown fields in y.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,tf],V. 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ϕ(a(z),z), we will identify thetotal derivative with the roman symbol d and thepartial derivative with the symbol ∂. It thus holds (12)ddzϕaz,z=∂ϕ∂adadz+∂ϕ∂zThis 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 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(13)dUdtP=QuP+TuP,where Qu is the power due to heat transfer and Tu 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(14a)Qu=∫PsqdA-∮∂Pq→·t→⊥dl(14b)Tu=∫PμLusLdA-∮∂PμRuh→R·t→⊥dlThe time variation of net internal energy U thus corresponds to the power expenditure of two external agents: a heat contribution Qu where sq is the heat supplied by external agents and q→ is the heat flux vector; a mass contribution Tu in which the scalar μβu denotes the change in specific energy provided by a unit supply ofmoles of species β=L,R.Figure 2
Notation.Since the geometry remains unchanged, one can define specific internal energyu per unit mass or per unit surface, since none of them changes during the process. We choose to define it per unit surface, namely,(15)UP=∫PudA.Standard application of the surface divergence theorem and of mass balances (2a)–(2c) leads from (14a) and (14b) to(16a)Qu=∫PsqdA-∫PdivΩq→dA(16b)Tu=∫PμLusLdA-∫PdivΩμRuh→RdAThe first law of thermodynamics is thus stated as follows:(17)∫PdudtdA=∫PsqdA-∫PdivΩq→dA-∫PdivΩμRuh→RdA+∫PμLusLdAIt must hold for any region P, since the latter is arbitrary. After simple algebra, the local form of the first principle thus reads(18)dudt=sq-divΩq→+μRu∂cR∂t+μLu∂cL∂t+μCu∂cC∂t-h→R·∇ΩμRu+μRu+μLu-μCuw1


#### 2.4.2. Entropy Balance Equations

The second law of thermodynamics represents the balance of the interplay among the internal entropy ofP and the entropy transferred in P due to mass exchange and heat transferred on P. The entropy balance for the problem at hand reads(19)dSdtP-dSirrdtP=QηP+TηP,where S is the net internal entropy of P, Sirr is the entropy produced inside P, Qη is the entropy per unit time due to heat transfer, and Tη is the entropy per unit time due to mass transfer. The individual contributions read(20a)Qη=∫PsqTdA-∮∂Pq→T·t→⊥dl,(20b)Tη=∫PμLηsLdA-∮∂PμRηh→R·t→⊥dl.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(21)dSirrdt≥0.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(22)∫Pddtη-sqT+divΩq→T-μLηsL+divΩμRηh→RdA≥0By multiplying per T≥0, replacing -sq+divΩq→ by means of the energy balance (18), and some simple algebra, the entropy imbalance becomes(23)∫PTdηdt-dudt-1Tq→·∇ΩT+∂cR∂tμRu-TμRη+∂cL∂tμLu-TμLη+∂cC∂tμCu-TμCη+Th→R·∇ΩμRη+μRu∂cR∂t+μLu∂cL∂t+μCu∂cC∂t-h→R·∇ΩμRu+μRu-TμRη+μLu-TμLη-μCu+TμCηw1dA≥0Denote with the symbol μβ the quantity(24)μβ=μβu-Tμβηand with the symbol A1 the following:(25)A1=-μR-μL+μCBy noting that(26)Th→R·∇ΩμRη=h→R·∇ΩTμRη-h→R·∇ΩTμRηone finally writes the entropy balance as(27)∫PTdηdt-dudt-1Tq→·∇ΩT+μR∂cR∂t+μL∂cL∂t+μC∂cC∂t-A1w1-h→R·∇ΩμR-h→R·∇ΩTμRηdA≥0


#### 2.4.3. Helmholtz Free Energy

The specific Helmholtz free energy is defined as(28)ψ=u-Tηand is taken as a function of temperature and concentrations, ψT,cR,cL,cC. It thus holds(29)Tdηdt-dudt=-dψdt-η∂T∂t=-∂ψ∂cL∂cL∂t-∂ψ∂cR∂cR∂t-∂ψ∂cC∂cC∂t-η+∂ψ∂T∂T∂twhich can be inserted in (27) to derive the entropy imbalance in the final form:(30)∫P-∂ψ∂cR+μR∂cR∂t+-∂ψ∂cL+μL∂cL∂t+-∂ψ∂cC+μC∂cC∂t-η+∂ψ∂T∂T∂tdA+∫P-1Tq_→·∇ΩT-A1w1-h→R·∇ΩμRdA≥0where q_→=q→+TμRηh→R


#### 2.4.4. Thermodynamic Restrictions

Inequality (30) must hold for any region P, since the latter was arbitrarily taken. Therefore, the following local inequality, usually termed after Clausius-Duhem, yields(31)-∂ψ∂cR+μR∂cR∂t+-∂ψ∂cL+μL∂cL∂t+-∂ψ∂cC+μC∂cC∂t-η+∂ψ∂T∂T∂t+-1Tq_→·∇ΩT-Aw1-h→R·∇ΩμR≥0This inequality must hold for any value of the time derivative of the temperature and of the concentrationscR, cL, and cC. 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(32)μR=∂ψ∂cR,μL=∂ψ∂cL,μC=∂ψ∂cC,η=-∂ψ∂TIn view of formula (32), the amount μβ declared in (24) acquires the meaning ofchemical potential and hence the term A(1) in (25) turns out to be theaffinity of the reaction (1).Equation (32) yields to the so-called Clausius-Plank inequality:(33)-1Tq_→·∇ΩT-A1w1-h→R·∇ΩμR≥0that splits under the assumptions of Curie’s principle and thermal equilibrium in the following set of inequalities:(34a)h→R·∇ΩμR≤0,(34b)A1w1≤0.


### 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:(35)ψcR,cL,cC=ψRcR+ψLcL+ψCcCThe 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 β=R,L,C:(36)ψβcβ=μβ0cβ+RTcβmaxϑβln⁡ϑβ+1-ϑβln⁡1-ϑβwhere ϑβ=cβ/cβmax is the ratio between the concentration and the saturation limit for each species. The chemical potential μβ can be written accordingly to (32) as(37)μβ=∂ψ∂cβ=μβ0+RTln⁡ϑβ-ln⁡1-ϑβA strategy to satisfy the thermodynamic restriction (34a) is to model the flux of receptors by Fickian-diffusion, which linearly correlates h→R to the gradient of its chemical potential μR:(38)h→R=-MRcR∇ΩμRby means of a positive definite mobility tensor MR. The following isotropic nonlinear [15] specialization for the mobility tensor MR(39)MRcR=uRcRmaxϑR1-ϑR1accounts for saturation. In formula (39), ϑR=cR/cRmax; cRmax is the saturation limit for receptors. ThemobilityuR>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 (cR=0) and the saturated (cR=cRmax) phases have vanishing mobilities. Neither the mobility uR nor the saturation concentration cRmax 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 (40)∇ΩμR=RT1cRmax1ϑR1-ϑR∇ΩcRFick’s Law (38) specializes as follows:(41)h→R=-DR∇ΩcRwhere DR=uRRT is the receptordiffusivity.


### 2.6. Chemical Kinetics

The chemical kinetics of reaction (1) is modeled via the law of mass action:(42)w1=kfϑL1-ϑLϑR1-ϑR-kbϑC1-ϑCAt chemical equilibrium, as w(1)=0 and A(1)=0, the concentrations obey the relation(43)kfkb=ϑCeq1-ϑCeq1-ϑReqϑReq1-ϑLeqϑLeq=exp⁡-μC0-μL0-μR0RT=Keq1which defines the constant of equilibrium Keq1 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 cC is related to the others by the equation A1=0, i.e., from (25) and (37)(44)A1=ΔG0+RTln⁡ϑC-ln⁡ϑR-ln⁡ϑL-ln⁡1-ϑC-ln⁡1-ϑL-ln⁡1-ϑR=0where ΔG0=μC0-μL0-μR0 is the standard Gibbs free energy. Far from saturation, when cβ≪cβmax,(45)cC=cRcLαhaving denoted with α the following constant:(46)α=cRmaxcLmaxcCmaxexp⁡ΔG0RT


### 2.7. Governing Equations

Making use of (45) and (7a), (7c) finally becomes(47)cLα+1∂cR∂t+cRα∂cL∂t+divΩ⁡-DR∇ΩcR=0By subtracting (7b) from (7a), it comes out that(48)∂cR∂t+divΩ⁡-DR∇ΩcR-∂cL∂t+sLx,t=0with sL from (5). These are the governing equations for cL, cR, whereas cC can be recovered from (45). Initial conditions read(49a)cRt=0=cR0(49b)cLt=0=0(49c)cCt=0=0


### 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 c^R a test function, reads(50)∫ΩcLα+1∂cR∂t+cRα∂cL∂t+divΩ⁡-DR∇ΩcRc^RdA=0Applying the divergence theorem over Ω, the former equation transforms as follows:(51a)∫ΩcLα+1∂cR∂tc^RdA+∫ΩcRα∂cL∂tc^RdA+DR∫Ω∇ΩcR·∇Ωc^RdA=0The weak form of (48), defining with c^L a test function, reads after easy algebra: (51b)∫Ωc^L∂cR∂tdA-∫Ωc^L∂cL∂tdA+DR∫Ω∇ΩcR·∇Ωc^LdA+∫Ωc^LsLx,tdA=0The 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 φi(x) and shape functions φj(x) and nodal unknowns that depend solely on time. The usual Einstein summation convention is taken henceforth for repeated indexes. (52)cR=φkRxcRkt,∇ΩcR=∇ΩφkRxcRkcL=φjLxcLjt,∇ΩcL=∇ΩφjLxcLjc^R=φiRx,∇Ωc^R=∇ΩφiRxcL^=φiLx,∇ΩcL^=∇ΩφiLxThe nonlinear ODEs read(53a)∫ΩφiRxφkRxφjLxdAcLjtα∂cRkt∂tdA+∫ΩφiRxφkRxdA∂cRkt∂t+∫ΩφiRxφkRxφjLxdAcRktα∂cLjt∂t+DR∫Ω∇ΩφkRx·∇ΩφiRxdAcRkt=0(53b)∫ΩφiLxφkRxdA∂cRkt∂t-∫ΩφiLxφjLxdA∂cLjt∂t+∫ΩφiLsLx,tdA+DR∫Ω∇ΩφkRx·∇ΩφiLxdAcRkt=0For 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(54)∂cRk∂t≃cRkt+Δt-cRktΔt,∂cLj∂t≃cLjt+Δt-cLjtΔt.We make recourse to the Backward Euler method that leads to the following nonlinear problem in cRk(t+Δt) and cLk(t+Δt):(55a)∫ΩφiRxφkRxφjLxdAcLjt+ΔtαcRkt+ΔtΔtdA-∫ΩφiRxφkRxφjLxdAcLjt+ΔtαcRktΔtdA+∫ΩφiRxφkRxdAcRkt+ΔtΔt-∫ΩφiRxφkRxdAcRktΔt+∫ΩφiRxφkRxφjLxdAcRkt+ΔtαcLjt+ΔtΔt-∫ΩφiRxφkRxφjLxdAcRkt+ΔtαcLjtΔt+DR∫Ω∇ΩφkRx·∇ΩφiRxdAcRkt+Δt=0(55b)∫ΩφiLxφkRxdAcRkt+ΔtΔt-∫ΩφiLxφkRxdAcRktΔt-∫ΩφiLxφjLxdAcLjt+ΔtΔt+∫ΩφiLxφjLxdAcLjtΔt+DR∫Ω∇ΩφkRx∇ΩφiLxdAcRkt+Δt+∫ΩφiLsLx,tdA=0


## 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 (R) and ligands (L) is described as a chemical reaction, which produces a receptor-ligand complex (C)(1)R+L⇄kbkfC,where kf and kb are the kinetic constants of the forward and backward reaction, respectively. The reaction rate w(1), measured in [mol/m3s], 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:(2a)∂cR∂t+divΩ⁡h→R+w1=sR,(2b)∂cL∂t+divΩ⁡h→L+w1=sL,(2c)∂cC∂t+divΩ⁡h→C-w1=sC.Symbols in (2a)–(2c) have the following meaning (concentrations cβ are defined in space and time, i.e., cβ=cβx→,t. The same holds for h→β, w1, and sβ. Functional dependence is specified when necessary only, to favor readability): cβ (with β=R,L,C ) is the molarity (i.e., the number of moles per unit area) of a generic species β; h→β 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; sβ is the rate in number of molecules per unit volume per unit time at which species β is generated by sources, and t 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.,(3)h→L=h→C=0→.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 n→ the cell membrane unit normal,(4)divΩ⁡h→R=n→·curln→×h→R=div⁡h→R-∇h→Rn→·n→.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 sL whose expression is calibrated from experimental evidence [14]. The following expression for sL in (2b) is taken [7]:(5)sLx,t=cL¯t¯Ht-xvHt¯-t+xvThe 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.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 functionsL in (5).In (5), H[-] is the Heaviside step function, cL¯=72 ligands/μm2 is the concentration of substrate-immobilized ligand available for reaction (1), tf is the time required for the complete mechanical deformation of the cell, v=πl/2tf is the velocity of mechanical deformation (assumed to be constant until tf), l is the cell radius, t¯≪tf is a parameter that identifies a finite time required for binding, x is the curvilinear abscissa of our simplified geometry, and t is the generic time. In view of (5), the supply of ligands at point x on the membrane remains zero until t<x/v; then, in the time span between t=x/v and t=x/v+t¯, it increases rapidly from zero to cL¯. We assume(6)sR=sC=0,since complex is provided by w(1) only, and receptors are not generated.In view of the above, mass balance equations (2a)–(2c) finally become(7a)∂cR∂t+divΩ⁡h→R+w1=0(7b)∂cL∂t+w1-sLx,t=0(7c)∂cC∂t-w1=0


## 2.3. Weak Form and Boundary Conditions

The weak formulation of balance equations (7a)–(7c) comes out after multiplication by a suitable set oftime 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:(8)∫Ωc^R∂cR∂t+divΩ⁡h→R+w1dA=∫Ωc^R∂cR∂tdA+∫ΩdivΩc^Rh→R-∇Ωc^R·h→RdA+∫Ωc^Rw1dA=∫Ωc^R∂cR∂tdA-∫Ω∇Ωc^R·h→RdA+∫Ωc^Rw1dA=0.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(9)∇Ωc^R=∇c^R-n→·∇c^Rn→with n→ 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 interval0,tf as(10)Findy∈V0,tfsuch that∂∂tby^,yt+ay^,yt=fy^∀y^∈Vwhere(11)by^,y=∫Ωc^RcR+c^LcL+c^CcCdA,ay^,y=-∫Ω∇Ωc^R·h→RdA+∫Ωc^R+c^L-c^Cw1dA,fy^=∫Ωc^LsLdAwith y=cR,cL,cC and y^=c^R,c^L,c^C. Column y collects the time-dependent unknown fields. Column y^ collects the steady-state test functions that correspond to the unknown fields in y.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,tf],V. 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ϕ(a(z),z), we will identify thetotal derivative with the roman symbol d and thepartial derivative with the symbol ∂. It thus holds (12)ddzϕaz,z=∂ϕ∂adadz+∂ϕ∂zThis 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 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(13)dUdtP=QuP+TuP,where Qu is the power due to heat transfer and Tu 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(14a)Qu=∫PsqdA-∮∂Pq→·t→⊥dl(14b)Tu=∫PμLusLdA-∮∂PμRuh→R·t→⊥dlThe time variation of net internal energy U thus corresponds to the power expenditure of two external agents: a heat contribution Qu where sq is the heat supplied by external agents and q→ is the heat flux vector; a mass contribution Tu in which the scalar μβu denotes the change in specific energy provided by a unit supply ofmoles of species β=L,R.Figure 2
Notation.Since the geometry remains unchanged, one can define specific internal energyu per unit mass or per unit surface, since none of them changes during the process. We choose to define it per unit surface, namely,(15)UP=∫PudA.Standard application of the surface divergence theorem and of mass balances (2a)–(2c) leads from (14a) and (14b) to(16a)Qu=∫PsqdA-∫PdivΩq→dA(16b)Tu=∫PμLusLdA-∫PdivΩμRuh→RdAThe first law of thermodynamics is thus stated as follows:(17)∫PdudtdA=∫PsqdA-∫PdivΩq→dA-∫PdivΩμRuh→RdA+∫PμLusLdAIt must hold for any region P, since the latter is arbitrary. After simple algebra, the local form of the first principle thus reads(18)dudt=sq-divΩq→+μRu∂cR∂t+μLu∂cL∂t+μCu∂cC∂t-h→R·∇ΩμRu+μRu+μLu-μCuw1


### 2.4.2. Entropy Balance Equations

The second law of thermodynamics represents the balance of the interplay among the internal entropy ofP and the entropy transferred in P due to mass exchange and heat transferred on P. The entropy balance for the problem at hand reads(19)dSdtP-dSirrdtP=QηP+TηP,where S is the net internal entropy of P, Sirr is the entropy produced inside P, Qη is the entropy per unit time due to heat transfer, and Tη is the entropy per unit time due to mass transfer. The individual contributions read(20a)Qη=∫PsqTdA-∮∂Pq→T·t→⊥dl,(20b)Tη=∫PμLηsLdA-∮∂PμRηh→R·t→⊥dl.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(21)dSirrdt≥0.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(22)∫Pddtη-sqT+divΩq→T-μLηsL+divΩμRηh→RdA≥0By multiplying per T≥0, replacing -sq+divΩq→ by means of the energy balance (18), and some simple algebra, the entropy imbalance becomes(23)∫PTdηdt-dudt-1Tq→·∇ΩT+∂cR∂tμRu-TμRη+∂cL∂tμLu-TμLη+∂cC∂tμCu-TμCη+Th→R·∇ΩμRη+μRu∂cR∂t+μLu∂cL∂t+μCu∂cC∂t-h→R·∇ΩμRu+μRu-TμRη+μLu-TμLη-μCu+TμCηw1dA≥0Denote with the symbol μβ the quantity(24)μβ=μβu-Tμβηand with the symbol A1 the following:(25)A1=-μR-μL+μCBy noting that(26)Th→R·∇ΩμRη=h→R·∇ΩTμRη-h→R·∇ΩTμRηone finally writes the entropy balance as(27)∫PTdηdt-dudt-1Tq→·∇ΩT+μR∂cR∂t+μL∂cL∂t+μC∂cC∂t-A1w1-h→R·∇ΩμR-h→R·∇ΩTμRηdA≥0


### 2.4.3. Helmholtz Free Energy

The specific Helmholtz free energy is defined as(28)ψ=u-Tηand is taken as a function of temperature and concentrations, ψT,cR,cL,cC. It thus holds(29)Tdηdt-dudt=-dψdt-η∂T∂t=-∂ψ∂cL∂cL∂t-∂ψ∂cR∂cR∂t-∂ψ∂cC∂cC∂t-η+∂ψ∂T∂T∂twhich can be inserted in (27) to derive the entropy imbalance in the final form:(30)∫P-∂ψ∂cR+μR∂cR∂t+-∂ψ∂cL+μL∂cL∂t+-∂ψ∂cC+μC∂cC∂t-η+∂ψ∂T∂T∂tdA+∫P-1Tq_→·∇ΩT-A1w1-h→R·∇ΩμRdA≥0where q_→=q→+TμRηh→R


### 2.4.4. Thermodynamic Restrictions

Inequality (30) must hold for any region P, since the latter was arbitrarily taken. Therefore, the following local inequality, usually termed after Clausius-Duhem, yields(31)-∂ψ∂cR+μR∂cR∂t+-∂ψ∂cL+μL∂cL∂t+-∂ψ∂cC+μC∂cC∂t-η+∂ψ∂T∂T∂t+-1Tq_→·∇ΩT-Aw1-h→R·∇ΩμR≥0This inequality must hold for any value of the time derivative of the temperature and of the concentrationscR, cL, and cC. 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(32)μR=∂ψ∂cR,μL=∂ψ∂cL,μC=∂ψ∂cC,η=-∂ψ∂TIn view of formula (32), the amount μβ declared in (24) acquires the meaning ofchemical potential and hence the term A(1) in (25) turns out to be theaffinity of the reaction (1).Equation (32) yields to the so-called Clausius-Plank inequality:(33)-1Tq_→·∇ΩT-A1w1-h→R·∇ΩμR≥0that splits under the assumptions of Curie’s principle and thermal equilibrium in the following set of inequalities:(34a)h→R·∇ΩμR≤0,(34b)A1w1≤0.


## 2.4.1. 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(13)dUdtP=QuP+TuP,where Qu is the power due to heat transfer and Tu 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(14a)Qu=∫PsqdA-∮∂Pq→·t→⊥dl(14b)Tu=∫PμLusLdA-∮∂PμRuh→R·t→⊥dlThe time variation of net internal energy U thus corresponds to the power expenditure of two external agents: a heat contribution Qu where sq is the heat supplied by external agents and q→ is the heat flux vector; a mass contribution Tu in which the scalar μβu denotes the change in specific energy provided by a unit supply ofmoles of species β=L,R.Figure 2
Notation.Since the geometry remains unchanged, one can define specific internal energyu per unit mass or per unit surface, since none of them changes during the process. We choose to define it per unit surface, namely,(15)UP=∫PudA.Standard application of the surface divergence theorem and of mass balances (2a)–(2c) leads from (14a) and (14b) to(16a)Qu=∫PsqdA-∫PdivΩq→dA(16b)Tu=∫PμLusLdA-∫PdivΩμRuh→RdAThe first law of thermodynamics is thus stated as follows:(17)∫PdudtdA=∫PsqdA-∫PdivΩq→dA-∫PdivΩμRuh→RdA+∫PμLusLdAIt must hold for any region P, since the latter is arbitrary. After simple algebra, the local form of the first principle thus reads(18)dudt=sq-divΩq→+μRu∂cR∂t+μLu∂cL∂t+μCu∂cC∂t-h→R·∇ΩμRu+μRu+μLu-μCuw1


## 2.4.2. Entropy Balance Equations

The second law of thermodynamics represents the balance of the interplay among the internal entropy ofP and the entropy transferred in P due to mass exchange and heat transferred on P. The entropy balance for the problem at hand reads(19)dSdtP-dSirrdtP=QηP+TηP,where S is the net internal entropy of P, Sirr is the entropy produced inside P, Qη is the entropy per unit time due to heat transfer, and Tη is the entropy per unit time due to mass transfer. The individual contributions read(20a)Qη=∫PsqTdA-∮∂Pq→T·t→⊥dl,(20b)Tη=∫PμLηsLdA-∮∂PμRηh→R·t→⊥dl.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(21)dSirrdt≥0.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(22)∫Pddtη-sqT+divΩq→T-μLηsL+divΩμRηh→RdA≥0By multiplying per T≥0, replacing -sq+divΩq→ by means of the energy balance (18), and some simple algebra, the entropy imbalance becomes(23)∫PTdηdt-dudt-1Tq→·∇ΩT+∂cR∂tμRu-TμRη+∂cL∂tμLu-TμLη+∂cC∂tμCu-TμCη+Th→R·∇ΩμRη+μRu∂cR∂t+μLu∂cL∂t+μCu∂cC∂t-h→R·∇ΩμRu+μRu-TμRη+μLu-TμLη-μCu+TμCηw1dA≥0Denote with the symbol μβ the quantity(24)μβ=μβu-Tμβηand with the symbol A1 the following:(25)A1=-μR-μL+μCBy noting that(26)Th→R·∇ΩμRη=h→R·∇ΩTμRη-h→R·∇ΩTμRηone finally writes the entropy balance as(27)∫PTdηdt-dudt-1Tq→·∇ΩT+μR∂cR∂t+μL∂cL∂t+μC∂cC∂t-A1w1-h→R·∇ΩμR-h→R·∇ΩTμRηdA≥0


## 2.4.3. Helmholtz Free Energy

The specific Helmholtz free energy is defined as(28)ψ=u-Tηand is taken as a function of temperature and concentrations, ψT,cR,cL,cC. It thus holds(29)Tdηdt-dudt=-dψdt-η∂T∂t=-∂ψ∂cL∂cL∂t-∂ψ∂cR∂cR∂t-∂ψ∂cC∂cC∂t-η+∂ψ∂T∂T∂twhich can be inserted in (27) to derive the entropy imbalance in the final form:(30)∫P-∂ψ∂cR+μR∂cR∂t+-∂ψ∂cL+μL∂cL∂t+-∂ψ∂cC+μC∂cC∂t-η+∂ψ∂T∂T∂tdA+∫P-1Tq_→·∇ΩT-A1w1-h→R·∇ΩμRdA≥0where q_→=q→+TμRηh→R


## 2.4.4. Thermodynamic Restrictions

Inequality (30) must hold for any region P, since the latter was arbitrarily taken. Therefore, the following local inequality, usually termed after Clausius-Duhem, yields(31)-∂ψ∂cR+μR∂cR∂t+-∂ψ∂cL+μL∂cL∂t+-∂ψ∂cC+μC∂cC∂t-η+∂ψ∂T∂T∂t+-1Tq_→·∇ΩT-Aw1-h→R·∇ΩμR≥0This inequality must hold for any value of the time derivative of the temperature and of the concentrationscR, cL, and cC. 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(32)μR=∂ψ∂cR,μL=∂ψ∂cL,μC=∂ψ∂cC,η=-∂ψ∂TIn view of formula (32), the amount μβ declared in (24) acquires the meaning ofchemical potential and hence the term A(1) in (25) turns out to be theaffinity of the reaction (1).Equation (32) yields to the so-called Clausius-Plank inequality:(33)-1Tq_→·∇ΩT-A1w1-h→R·∇ΩμR≥0that splits under the assumptions of Curie’s principle and thermal equilibrium in the following set of inequalities:(34a)h→R·∇ΩμR≤0,(34b)A1w1≤0.


## 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:(35)ψcR,cL,cC=ψRcR+ψLcL+ψCcCThe 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 β=R,L,C:(36)ψβcβ=μβ0cβ+RTcβmaxϑβln⁡ϑβ+1-ϑβln⁡1-ϑβwhere ϑβ=cβ/cβmax is the ratio between the concentration and the saturation limit for each species. The chemical potential μβ can be written accordingly to (32) as(37)μβ=∂ψ∂cβ=μβ0+RTln⁡ϑβ-ln⁡1-ϑβA strategy to satisfy the thermodynamic restriction (34a) is to model the flux of receptors by Fickian-diffusion, which linearly correlates h→R to the gradient of its chemical potential μR:(38)h→R=-MRcR∇ΩμRby means of a positive definite mobility tensor MR. The following isotropic nonlinear [15] specialization for the mobility tensor MR(39)MRcR=uRcRmaxϑR1-ϑR1accounts for saturation. In formula (39), ϑR=cR/cRmax; cRmax is the saturation limit for receptors. ThemobilityuR>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 (cR=0) and the saturated (cR=cRmax) phases have vanishing mobilities. Neither the mobility uR nor the saturation concentration cRmax 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 (40)∇ΩμR=RT1cRmax1ϑR1-ϑR∇ΩcRFick’s Law (38) specializes as follows:(41)h→R=-DR∇ΩcRwhere DR=uRRT is the receptordiffusivity.


## 2.6. Chemical Kinetics

The chemical kinetics of reaction (1) is modeled via the law of mass action:(42)w1=kfϑL1-ϑLϑR1-ϑR-kbϑC1-ϑCAt chemical equilibrium, as w(1)=0 and A(1)=0, the concentrations obey the relation(43)kfkb=ϑCeq1-ϑCeq1-ϑReqϑReq1-ϑLeqϑLeq=exp⁡-μC0-μL0-μR0RT=Keq1which defines the constant of equilibrium Keq1 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 cC is related to the others by the equation A1=0, i.e., from (25) and (37)(44)A1=ΔG0+RTln⁡ϑC-ln⁡ϑR-ln⁡ϑL-ln⁡1-ϑC-ln⁡1-ϑL-ln⁡1-ϑR=0where ΔG0=μC0-μL0-μR0 is the standard Gibbs free energy. Far from saturation, when cβ≪cβmax,(45)cC=cRcLαhaving denoted with α the following constant:(46)α=cRmaxcLmaxcCmaxexp⁡ΔG0RT


## 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 cC is related to the others by the equation A1=0, i.e., from (25) and (37)(44)A1=ΔG0+RTln⁡ϑC-ln⁡ϑR-ln⁡ϑL-ln⁡1-ϑC-ln⁡1-ϑL-ln⁡1-ϑR=0where ΔG0=μC0-μL0-μR0 is the standard Gibbs free energy. Far from saturation, when cβ≪cβmax,(45)cC=cRcLαhaving denoted with α the following constant:(46)α=cRmaxcLmaxcCmaxexp⁡ΔG0RT


## 2.7. Governing Equations

Making use of (45) and (7a), (7c) finally becomes(47)cLα+1∂cR∂t+cRα∂cL∂t+divΩ⁡-DR∇ΩcR=0By subtracting (7b) from (7a), it comes out that(48)∂cR∂t+divΩ⁡-DR∇ΩcR-∂cL∂t+sLx,t=0with sL from (5). These are the governing equations for cL, cR, whereas cC can be recovered from (45). Initial conditions read(49a)cRt=0=cR0(49b)cLt=0=0(49c)cCt=0=0


## 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 c^R a test function, reads(50)∫ΩcLα+1∂cR∂t+cRα∂cL∂t+divΩ⁡-DR∇ΩcRc^RdA=0Applying the divergence theorem over Ω, the former equation transforms as follows:(51a)∫ΩcLα+1∂cR∂tc^RdA+∫ΩcRα∂cL∂tc^RdA+DR∫Ω∇ΩcR·∇Ωc^RdA=0The weak form of (48), defining with c^L a test function, reads after easy algebra: (51b)∫Ωc^L∂cR∂tdA-∫Ωc^L∂cL∂tdA+DR∫Ω∇ΩcR·∇Ωc^LdA+∫Ωc^LsLx,tdA=0The 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 φi(x) and shape functions φj(x) and nodal unknowns that depend solely on time. The usual Einstein summation convention is taken henceforth for repeated indexes. (52)cR=φkRxcRkt,∇ΩcR=∇ΩφkRxcRkcL=φjLxcLjt,∇ΩcL=∇ΩφjLxcLjc^R=φiRx,∇Ωc^R=∇ΩφiRxcL^=φiLx,∇ΩcL^=∇ΩφiLxThe nonlinear ODEs read(53a)∫ΩφiRxφkRxφjLxdAcLjtα∂cRkt∂tdA+∫ΩφiRxφkRxdA∂cRkt∂t+∫ΩφiRxφkRxφjLxdAcRktα∂cLjt∂t+DR∫Ω∇ΩφkRx·∇ΩφiRxdAcRkt=0(53b)∫ΩφiLxφkRxdA∂cRkt∂t-∫ΩφiLxφjLxdA∂cLjt∂t+∫ΩφiLsLx,tdA+DR∫Ω∇ΩφkRx·∇ΩφiLxdAcRkt=0For 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(54)∂cRk∂t≃cRkt+Δt-cRktΔt,∂cLj∂t≃cLjt+Δt-cLjtΔt.We make recourse to the Backward Euler method that leads to the following nonlinear problem in cRk(t+Δt) and cLk(t+Δt):(55a)∫ΩφiRxφkRxφjLxdAcLjt+ΔtαcRkt+ΔtΔtdA-∫ΩφiRxφkRxφjLxdAcLjt+ΔtαcRktΔtdA+∫ΩφiRxφkRxdAcRkt+ΔtΔt-∫ΩφiRxφkRxdAcRktΔt+∫ΩφiRxφkRxφjLxdAcRkt+ΔtαcLjt+ΔtΔt-∫ΩφiRxφkRxφjLxdAcRkt+ΔtαcLjtΔt+DR∫Ω∇ΩφkRx·∇ΩφiRxdAcRkt+Δt=0(55b)∫ΩφiLxφkRxdAcRkt+ΔtΔt-∫ΩφiLxφkRxdAcRktΔt-∫ΩφiLxφjLxdAcLjt+ΔtΔt+∫ΩφiLxφjLxdAcLjtΔt+DR∫Ω∇ΩφkRx∇ΩφiLxdAcRkt+Δt+∫ΩφiLsLx,tdA=0


## 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 thedeal.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 tf = 600 s as the time for completion of the mechanical deformation and t¯=1s 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 l=20μm was calculated from the measure of radius of 50 endothelial cells (ECs); receptor diffusivity DR=0.198μm2/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 cR0=4.8molecules/μm2. The kinetic parameters ligand/receptor were measured by surface plasmon resonance (SPR) measurements and the equilibrium constant turned out to be Keq(1)=354059. Finally, it has been taken cLmax=16000molecules/μm2 and cRmax=cCmax.The initial discretization on the membrane is depicted in Figure3(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 thedeal.ii framework, we used it for the problem at hand. The mesh at time t=85s is depicted in Figure 3(b). It consists of 14058 quadrilaterals, graded around the surface of contact between the cell and the substrate.Figure 3
Initial (a) and locally refined (b) meshes.
(a) (b)Figure4 depicts the resulting evolution of the concentration of receptors in space and time. Figure 4(a) depicts a snapshot of cR at t=85s 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 cR0=4.8receptors/μm2. After 60 s, the concentration profile is perturbed at the bottom of the cell (s=0), and such a perturbation “propagates” with time. This event can be clearly seen in Figure 4(a) at t=85s. At t=600s, 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 sL 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/μm2 and a steady state has not yet been reached.Figure 4
Time evolution of the spatial concentrationcR of free VEGFR2 along the cell membrane. (a) Snapshot at t=85 s on a half of the membrane, with identification of points A, B, C, D, and E. (b) One-dimensional description of cR at different times. Each curve plots the distribution of free receptors at different times t=60n, with n = 0, 1, 2, …,120 s from the beginning of the experiment at t=0 to the final time t=7200 s. Points A, B, C, D, and E correspond to (a).
(a) (b)Immediately after the membrane gets in contact with the substrate, the evolution ofcR 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 t=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 Figure 5, showing that mass is conserved in our numerical validation.Figure 5
Time evolution of the overall amount of bound VEGFR2-gremlin 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 timet=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.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 Section2.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 (Section2.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.



---

*Source: 4705472-2018-08-26.xml*