A three-dimensional finite element model of cAMP signals

This paper presents a three-dimensional finite element model for cyclic adenosine monophosphate (cAMP) signaling. Governing equations for the synthesis, diffusion, and degradation of cAMP were numerically implemented using the finite element method. Simulated results were displayed as time course plots of cAMP concentrations at selected nodes within the discretized geometry. The validity of the finite element model was assessed by comparing simulated results against analytical or other numerical solutions of cAMP concentration distribution for a spherical cellular volume. An endothelial cell was also simulated using its discretized geometry obtained from microscopic cellular cross-sectional images. Simulated solutions using the spherical cellular volume produced near identical cAMP concentration plots to the analytical solutions and were in good agreements with numerical results obtained from VCell, an existing software package for modeling cell biological systems. The validated 3-D finite element model was then employed to simulate the cAMP signaling pathway within a pulmonary microvascular endothelial cell geometry.


Introduction
Although enzymatic control of the cyclic adenosine monophosphate (cAMP) signaling pathway is well understood, the way in which information is encoded within cAMP signals is not as well defined. cAMP is a second messenger in intracellular signaling involved in This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). the regulation of numerous cellular functions. These functions include cellular proliferation, differentiation, and gene expression [1,2]. Synthesis of cAMP occurs typically at the plasma membrane in response to stimulated adenylyl cyclase (AC) activity. Degradation of cAMP is due to phosphodiesterase (PDE) enzyme activity. There is evidence of correlations between cAMP levels and disease. For example, in vitro infections with HIV have shown T cells with higher concentrations of cAMP [3]. Pharmaceuticals regulating cAMP levels have been developed to treat diseases such as diabetes, hypertension, and asthma [3][4][5][6]. Phosphodiesterase inhibitors, in particular, have been proven as successful treatments for central nervous system damage [7].
With only rudimentary knowledge of information encoding within intracellular signals, we have limited understanding of the subcellular localization, kinetics, and frequency of cAMP signals [8,9]. Numerical models for predicting spatial distributions of cAMP and other second messengers have been developed. A fourth order Runge-Kutta MATLAB script was used to simulate cAMP distributions near the plasma membrane of HEK-293 cells [10][11][12]. A large-scale stochastic model was also developed to model the subcellular localization of cAMP signals in HEK-293 cells using the software NeuroRD [13]. Both the models and study outcomes have been summarized [14]. The program Virtual Cell [15], or VCell, is an intracellular modeling software developed by the University of Connecticut. It uses the Finite Volume Method to simulate deterministic partial differential equation systems onto uploaded cellular geometries [15].
This work focuses on the implementation of the finite element method (FEM) for threedimensional modeling of the cAMP signaling pathway. The finite element method is a numerical engineering technique useful in complex geometries. It relies on the discretization of the geometry and the conversion of the partial differential equations governing the system into a system of algebraic equations. A two-dimensional finite element model had been developed and reported in our previous work [16] where more information about related works on cAMP intracellular signaling can be found.
The three-dimensional finite element model used equations governing the synthesis, diffusion, and degradation of cAMP. Four node tetrahedral elements were used to discretize the cellular geometries. Simulations were then run with initial conditions and parameters to model the spatial and temporal dynamics of the cAMP signaling pathway. Both spherical and realistic cellular geometries have been modeled. Cellular geometries were obtained from endothelial cross-sectional cell image sequences that were converted into a solid, finite element meshed volume. The finite element model ran off MATLAB scripts.
The 3-D finite element simulations were validated through comparison of available analytical solutions and of resulting simulations previously done using the software VCell [17]. The validated 3-D FEM model was then utilized for simulating cAMP intracellular signaling within a pulmonary microvascular endothelial cell (PMVEC).
In a previous work from our group, a 3-D FEM model was developed [18]. However, that model was based upon a linear approximation of the Michaelis-Menten enzyme kinetics about the initial value of cAMP concentration. As a result, the model in [18] is only accurate for low cAMP concentrations. On the other hand, the 3-D FEM model proposed in this work utilized a quasi-linearization of the Michaelis-Menten kinetics. This technique is based on a linear approximation about the cAMP concentration level at any given time. Despite the fact that the quasi-linearization requires some iterations to converge, the technique produces an accurate estimate of PDE reactions at any levels of cAMP concentration. In addition, the 3-D model developed in this work is able to simulate a variety of AC and PDE activities anywhere within the cell, as well as the instants when these activities start.

Governing equation
For 3-D cAMP intracellular signaling, its governing equation describing cAMP synthesis, diffusion and degradation is generally given as (see [19,20]) where C = C(t, x, y, z) is the cAMP concentration at time t and location (x,y,z), D is the diffusion coefficient, E AC is the cAMP synthesis function, M(C) is the Michaelis-Menten cAMP degradation function, and t s and t d are the time points when AC activity (cAMP synthesis) and PDE activity (cAMP degradation) start taking place, respectively. In this work, PDE activity is assumed to occur after or at the same time of AC activity t s ≤ t d .
Under the steady-state assumption, where V max is the maximum cAMP hydrolysis rate, and K M is the Michaelis-Menten constant for cAMP binding to PDE.
As M(C) is a nonlinear function in C, an iterative method will need to be employed in the FEA of the model given by Eq. (1). The iterative technique adopted in this work is based on a quasi-linearization of M(C) about C = C 1 : To accurately evaluate Eq. (2) using its quasi-linearization (3), at a given time and location, the iteration must be carried out until the difference between the predicted C 1 and the solution C resulting from solving the third equation of system (1) satisfies a chosen convergence criterion.
The possible boundary conditions are  (4) where n is the normal vector to the boundary surface and β is a constant. For the proposed FEA model, this concentration flux β can be used to simulate AC activities in the plasmalemmal or perinuclear region of a cell.

Finite element implementation of the governing equation
To use the Galerkin approximation, we first discretize the cellular geometry into a number of elements. For each element, we multiply Eq. (1) by the shape functions N i (i = 1, 2, …, n, where n is the number of nodes of the chosen type of element) selected as weighting functions and then integrate it over the volume V of the element as follows: where ∇ 2 is the Laplacian, and By using Gauss's divergence theorem on the diffusion term and the boundary conditions (see, e.g., [21] for more details), one obtains the following weak form in 3-D: Warren et al. where S n is the face of the element over which its concentration flux is specified.
In this equation, the concentration C of cAMP is interpolated over the element from the nodal values c 1 , c 2 , …, c n using the shape functions N 1 , N 2 , …, N n as follows: Thus the time derivative of C and its gradients are given by  (10) Substitution of Eqs. (8)-(10) into Eq. (7) results in, where (13) [κ] = D 1 0 0 0 1 0 0 0 1 (14) Warren et al. Page 5 Forces Mech. Author manuscript; available in PMC 2022 January 20.
In this work, the four-node tetrahedral element (see Fig. 1) was selected for a 3-D numerical implementation of Eq. (11). The shape function matrix of this element is known to be (19) where ξ, η and ζ are the coordinates of the natural coordinate system. By using these shape functions in Eqs. (12), (13), (15)-(18), one gets where B = x 2 − x 1 y 2 − y 1 z 2 − z 1 In the above equations, (x i , y i , z i ) are the nodal coordinates of the four-node tetrahedral element under consideration, V is the volume of the element, and A ijk is the area of face i-j-k of the element.
At this point, the matrices and vectors in Eqs. (20)- (26) for all elements need to be expanded to the structure/model size before they can be assembled to obtain the global version of the first-order differential Eq. (11).
By using the time integration method [22], the vector of unknown concentrations {c} i+1 at time t i+1 can be found from {c} i at time t i as where Δt = t i + 1 − t i is the time step and in this work, Galerkin's implicit method (γ = 2/3) was chosen as the method is known to be unconditionally stable (no restriction on Δt for obtaining a stable solution, [22]). If some nodal concentrations are prescribed, these boundary conditions must be applied to Eq. (28) to obtain a reduced system of linear equations that contains the vector {c r } i+1 of only unknown nodal concentrations at time t i+1 .

The equation system (27) is simply a linear system of algebraic equations of the form
For each step time of the time integration method described in Eq. (27), as mentioned before, an iterative process must be utilized for the quasi-linearization of the Michaelis-Menten model. To be specific, the concentration solution {c} i from the previous step time t i will be used as an initial guess for C 1 employed in evaluating [K 3 ] and {R 2 } for calculations at step time t i+1 (see Eqs. (6), (23) and (26)). The solution {c*} i+1 resulting from using these [K 3 ] and {R 2 } is expected to be a better guess for C 1 and this process should be repeated until a chosen convergence criteria is met.

Validation of the proposed FEA model
The 3-D FEA model developed was validated against some available analytical solutions and published data. The FEA simulations employed the same cAMP signaling data as in [17] for PMVECs, i.e., E AC = 0.1412 μM/s for synthesis, D = 0.3 to 300 μm 2 /s for diffusion, V max = 0.295 μM/s, K M = 2 μM for degradation, and C o = 0.05 μM for initial conditions.
In addition, no cAMP signal is assumed to be transported between the cytoplasm and the nucleus.
The validation was first conducted against some analytical solutions. For this purpose, FEA results were obtained using a spherical cell model. The radii of the cell and the nucleus are R o = 9.34 μm and R i = 5.26 μm, respectively. The validated 3-D FEA model was then applied to simulate intra cellular cAMP signaling within a cultured PMVEC model. Fig. 3(a) and 3(b) depict the tetrahedral meshes with 4,983 nodes and 29,069 elements for the PMVEC model.
The FEA meshes shown in Fig. 2(a), 2(b), 3(a) and 3(b) are those satisfying the mesh convergence tests presented in section 4.1.2. These meshes also have good quality measured by their aspect ratio and solid angle as depicted in Table 1 [24]. For each quality measure, its best and worst possible values are 1 and 0, respectively. The mesh quality for the spherical cell model is better than the PMVEC model as the geometry of the latter is much more complex.  Fig. 4(a) and 4(b) confirm the advantage of Galerkin's implicit method as the cAMP responses remain stable with different time steps selected, namely, Δt = 10, 5 and 1 sec. For the spherical cell model using mesh 3 (see Table  2 and section 4.1.2), the numerical results are mostly unchanged and accurate even with the use of a very coarse time step of 10 sec. For the PMVEC model with a complex geometry (934 nodes and 5,039 elements), the error of the numerical results is visible at large time steps. However, the accuracy of the FEA result in case Δt = 10 sec is still acceptable. As expected, the results converge to the analytical solution as Δt is reduced.

Mesh convergence tests-A series of meshes with increasing mesh density
were employed for the mesh convergence tests. Table 2 shows the information for the two coarsest and two finest meshes in the series for the spherical cell model.
For compartmental models (see sections 4.2.1 and 4.2.2) the outputs are almost independent of the mesh density. As a result, the use of mesh 1 (see Table 2 Fig. 5(a) demonstrates this mesh independence behavior of compartmental models. In this test, the coarse mesh has 934 nodes and 5,039 tetrahedral elements while the fine mesh is shown in Fig. 3(a) and 3(b).
For non-compartmental models, such as those presented in section 4.3, the outputs are sensitive to the mesh density. Fig. 5

Validation against analytical solutions
There is no general analytical solution to the governing Eq. (1). However, analytical solutions are available for some particular cases involving simple cellular geometries or uniform loading conditions. These particular cases are considered in this section to partially verify the proposed FEA model.

Verification of the implementation of the Michaelis-Menten equation-
If there is no concentration flux (β = 0) at the subplasmalemmal and perinuclear regions, and cAMP synthesis and/or degradation occur uniformly in the cytosolic region, then there is no diffusion of cAMP which means cAMP concentration at any cellular location is the same at a given time. This is known as a compartmental model for which the resulting concentration C is independent of cellular location (x, y, z), cellular geometries and the diffusion coefficient D which means that C is a function of time only. Hence, for t ≥ t d , the partial differential Eq. (1) becomes the following ordinary differential equation for finding C(t): The analytical solution for this equation is given by Eq. (33) in the Appendix.
In this section, Eq. (29) was used to validate the FEA implementation of the iterative process for the quasi-linearization of the Michaelis-Menten Eq. (2). To this end, a compartmental model for which no AC activity (E AC = 0) is present and PDE activity is uniformly distributed within the cytosol from t d = 0 was studied. The concentration responses for four different scenarios of PDE activities are depicted in Fig. 6(a) and 6(b) where the FEA results agree perfectly with the analytical solutions given by Eq. (35). The outputs confirm the compartmental behavior of the FEA model as they are independent of the cellular position, the cellular geometry (PMVEC or spherical cell) and the diffusion coefficient used. As predicted, an increase of V max with respect to K M (Fig. 6(a)) speeds up the cAMP degradation, while an elevation of K M with respect to V max slows down the degradation.
As with the 2-D FEA model [16], the following biological behavior can be observed for compartmental models studied in this section: cAMP concentration responses are the same as long as the ratio of V max /K M is the same. This can be confirmed by noticing that the curve for V max = 0.118 μM/s and K M = 2 μM in Fig. 6(a) and the curve for V max = 0.295 μM/s and K M = 5 μM in Fig. 6(b) are identical as they have the same ratio V max /K M = 0.059.

Verification of the implementation of cAMP synthesis-Eq.
(29) was also employed to validate the implementation of cAMP synthesis in the proposed FEA model. This was done by using a compartmental scenario where the plasma membrane and the perinuclear region are under flux-free boundary conditions, PDE activity is absent (V max = 0) and AC activity is uniformly occurred inside the cytosol from t s = 0. The analytical solution for the compartmental output for this case is given by Eq. (36) Three AC activities characterized by three synthesis rates E AC = 0.11, 0.1412 and 0.17 μM/s were investigated. Very good agreement between the FEA results and the analytical solutions can be seen in Fig. 7(a).
In a next step, a further validation involving a compartmental case where both AC and PDE activities take place at the same time was considered. This was done by adding to the three AC activities a uniform PDE activity (V max = 0.295 μM/s, K M = 2 μM) within the cytosol from t d = 0. Fig. 7(b) shows an excellent agreement between the FEA results and the analytical counterparts obtained by employing Eq. (33). As expected, a lower amount of cAMP synthesis rate results in a shorter time before the steady state of concentrations being reached.

Verification of the implementation of cAMP diffusion-To validate the
proposed FEA model in terms of numerical implementation of cAMP diffusion, two simulations involving different boundary conditions at the subplasmalemmal and perinuclear regions were studied. Both simulations were conducted on the spherical cell geometry in Fig. 2(a) subjected to a diffusion coefficient D = 30 μm 2 /s and an absence of AC and PDE activities within the cytosol (E AC = 0 and V max = 0).
For the first simulation, cAMP concentrations of 4 μM and 1 μM were prescribed over the plasma membrane and perinuclear region, respectively. The FEA results for the time response of cAMP concentration on two spherical surfaces of radii R = 7.468 μm and R = 7.112 μm within the cytosol are favorably compared against the steady-state analytical results as shown in Fig. 8(a).
For the second simulation, concentration boundary conditions (C = 3 μM or C = 5 μM) were applied at the plasma membrane while flux-free boundary conditions were imposed over the perinuclear region. Because of these boundary conditions, the cAMP concentration at any locations in the subplasmalemmal region will remain at the applied value (3 μM or 5 μM) while the concentration at other locations in the cytosolic and perinuclear regions will rise from the initial value of 0.05 μM to the steady-state value of 3 μM or 5 μM. Due to this behavior, it is sufficient to show the FEA result for a representative location within the cytosolic region.

Validation against other numerical technique
In this section, two simulations involving non-compartmental models for the spherical cell geometry (Fig. 2(a) and 2(b)) were used to verify the proposed FEA technique as a whole by comparing the FEA results with those obtained from the finite volume technique implemented within the VCell software [15].
In both simulations, AC activity was uniformly distributed at the plasma membrane and PDE activity was uniformly distributed in the cytosol. The first simulation employed a diffusion coefficient of 3 μm 2 /s while that value chosen for the second simulation was 0.3 μm 2 /s. Note that these simulations were previously run using the VCell software and the results were reported in [17].
The uniform AC distribution at the subplasmalemmal region was modeled using a positive flux on the outer surface of the spherical cell geometry. As in [17], β was determined under the assumption of the same total AC activity produced either over the plasma membrane or inside the cytosolic region, i.e., where E AC = 0.1412 μM/s, V c and S p are the volume of the cytosolic region and surface area of the plasma membrane of the spherical cell geometry shown in Fig. 2(a).
Therefore, Fig. 9(a) and 9(b) show a comparison between the FEA and VCell results for the first (D = 3 μm 2 /s) and second (D = 0.3 μm 2 /s) simulations, respectively. In each figure, the time responses of cAMP concentration on the plasma membrane (R o ), on the spherical surface of radius R = 7.8μm and at the perinuclear region (R i ) are plotted. There is reasonable agreement between the two results (FEA vs VCell), but not perfect. For the second simulation (see Fig. 9(b)), some noticeable discrepancy occurs at the perinuclear region in the early time of the process (t < 10 s). However, the FEA curve during this early time period makes sense as the perinuclear region is far from the source of synthesis while D is small (cAMP spread is slow): cAMP concentration at the perinuclear region should first drop from the initial value of 0.05 μM (due to PDE activity) before AC activity reaches this region and raises C to its steady-state value. The absence of this concentration drop on the VCell curve may be explained by the fact that the VCell simulation in [17] used a larger time step (Δt = 10 s vs Δt = 2 s employed in the FEA simulation) which prevented it to capture the concentration drop within the first 10 seconds of the process.

FEA simulations on the 3-D cultured PMVEC geometry
The proposed FEA model previously validated was used in this section to simulate cAMP cellular signaling within the complex 3-D cultured PMVEC geometry shown in Fig. 3(a) and 3(b). Two simulations with and without delayed start of AC and PDE activities were considered herein. The time responses of cAMP concentration at three typical cellular locations were sought. These locations are represented by nodes 2438, 2460 and 2437 at the subplasmalemmal, cytosolic and perinuclear regions, respectively. These nodes have the same z-coordinate of 4.6 μm. The location of these nodes in the cross section z = 4.6 μm is depicted in Fig. 10.

Uniform AC and PDE activities-
In this simulation, AC activity uniformly generated from the plasma membrane was modeled using a positive flux β = 2 μM·μm/s. PDE activity (V max = 0.295 μM/s, K M = 2 μM) was uniformly taken place throughout the cytosol. A diffusion coefficient D = 10 μm 2 /s was selected. There was no delayed start of AC and PDE activities, t s = t d = 0s. Fig. 11 shows the time responses of cAMP concentration at three chosen cellular locations. As expected, all the three curves start from the initial conditions and then approach their respective steady-state values. The closer the node to the source of synthesis, the higher the amount of steady-state concentration. As nodes 2460 and 2437 are at some distances from the plasma membrane where AC activity occurs, their curves exhibit an initial drop of cAMP concentration due to cAMP degradation before they recover and reach steady-state level.

Uniform AC and PDE activities with delayed start-
The only difference between this and the previous simulation is the delayed start of cAMP synthesis and degradation: AC and PDE activities were assumed to start at t s = 10 s and t d = 20 s, respectively. As a result, the concentrations at the three cellular locations remained at the initial level of 0.05 μM until t = 10 s when they started to rise (see Fig. 12). Due to cAMP synthesis without degradation between t = 10 and 20 s, there is a linear portion within this time period on each of the three curves depicted in Fig. 12. As soon as PDE activity started at t = 20 s, the slopes of the three curves decreased and the concentration started approaching their respective steady-state values.

Conclusion
A developed 3-D finite element model for cAMP intracellular signaling was presented in this paper. The model produced time course plots of cAMP concentrations at selected nodes using equations governing the synthesis, diffusion, and degradation of the second messenger. The finite element model is capable of simulating multiple AC and PDE activity scenarios. It can also confine cAMP synthesis and degradation to certain areas of the cell such as uniformly distributing both AC and PDE activity within the cytosol and/or bounding AC activity to the plasma membrane with a boundary flux condition. Initial cAMP concentrations, diffusion coefficients, and simulation times are also modifiable variables in the model.
The time course cAMP concentration plots simulated with the spherical cell geometry produced almost identical curves to the analytical literature solutions. Due to the use of two different numerical techniques (FEM vs FVM), there was a minor but uniform discrepancy between the simulated FEA and VCell results when a boundary flux β was specified to contain AC activity to the plasma membrane while PDE activity was uniformly distributed within the cytosol. Finally, the validated 3-D FEM model was successfully used to simulate cAMP intracellular signaling within a cultured PMVEC cell subjected to different distributions of AC and PDE activities.
In general, similar to a conclusion in our previous 2-D work [16], 3-D simulations also showed that sustained cAMP gradients can be formed within endothelial cells which is consistent with those observed in rat PMVECs [26]. These data demonstrate that the finite element modeling approach is a viable approach for modeling second messenger signaling systems in four dimensions (x, y, z, t). While the simulation data employed in this work were acquired from PMVECs, applications of the proposed FEA model to other cell types are possible.
where E AC , V max , K M and C 0 are constants, is given by where W() is the Lambert W function, and If E AC = 0, the solution reduces to If V max = 0, the solution is simply given by

Fig. 1.
A four-node tetrahedral element.         FEA solution for the time responses of cAMP concentration at nodes 2438 (on the plasma membrane), 2460 (in the cytosol) and 2437 (at the perinuclear region) in response to AC activity (β = 2 μM·μm/s) uniformly distributed over the plasma membrane at t s = 0, PDE activity uniformly distributed in the cytosol at t d = 0, and D = 10 μm 2 /s.  FEA solution for the time responses of cAMP concentration at nodes 2438 (on the plasma membrane), 2460 (in the cytosol) and 2437 (at the perinuclear region) in response to AC activity (β = 2 μM·μm/s) uniformly distributed along the plasma membrane at t s = 10 s, PDE activity uniformly distributed in the cytosol at t d = 20 s, and D = 10 μm 2 /s.