Computational Model of Ca2+ Wave Propagation in Human Retinal Pigment Epithelial ARPE-19 Cells

Objective Computational models of calcium (Ca2+) signaling have been constructed for several cell types. There are, however, no such models for retinal pigment epithelium (RPE). Our aim was to construct a Ca2+ signaling model for RPE based on our experimental data of mechanically induced Ca2+ wave in the in vitro model of RPE, the ARPE-19 monolayer. Methods We combined six essential Ca2+ signaling components into a model: stretch-sensitive Ca2+ channels (SSCCs), P2Y2 receptors, IP3 receptors, ryanodine receptors, Ca2+ pumps, and gap junctions. The cells in our epithelial model are connected to each other to enable transport of signaling molecules. Parameterization was done by tuning the above model components so that the simulated Ca2+ waves reproduced our control experimental data and data where gap junctions were blocked. Results Our model was able to explain Ca2+ signaling in ARPE-19 cells, and the basic mechanism was found to be as follows: 1) Cells near the stimulus site are likely to conduct Ca2+ through plasma membrane SSCCs and gap junctions conduct the Ca2+ and IP3 between cells further away. 2) Most likely the stimulated cell secretes ligand to the extracellular space where the ligand diffusion mediates the Ca2+ signal so that the ligand concentration decreases with distance. 3) The phosphorylation of the IP3 receptor defines the cell’s sensitivity to the extracellular ligand attenuating the Ca2+ signal in the distance. Conclusions The developed model was able to simulate an array of experimental data including drug effects. Furthermore, our simulations predict that suramin may interfere ligand binding on P2Y2 receptors or accelerate P2Y2 receptor phosphorylation, which may partially be the reason for Ca2+ wave attenuation by suramin. Being the first RPE Ca2+ signaling model created based on experimental data on ARPE-19 cell line, the model offers a platform for further modeling of native RPE functions.


Introduction
Epithelial tissue covers and lines all internal and external body surfaces. These cell layers have multiple functions depending on their location, and many of these functions are controlled by Ca 2+ activity [1]. Retinal pigment epithelium (RPE), a monolayer of pigmented polarized cells, is crucial for the maintenance of visual functions. Located in the back of the eye between photoreceptors and choriocapillaries, RPE forms a vital part of the blood-retinal barrier (BRB) [2]. The physiology of RPE is tightly coupled with the activity of the various ion channels, such as Ca 2+ channels that are associated with several important RPE functions including transepithelial transport of ions and water, dark adaption of photoreceptor activity, phagocytosis, secretion, and differentiation [3]. In RPE, as well as in other epithelia, local deformation of the cell membrane initiates a significant Ca 2+ wave [4][5][6]. Such deformation of the cell membrane can occur in clinically important pathological conditions such as retinal tear resulting from complications after photodynamic therapy [7], intravitreal bevacizumab injection [8], or intravitreal pegaptanib injection [9]. Intercellular Ca 2+ signaling is also linked to the initial stages of wound repair: excessive mechanical stimulation causes cell death and thus initiates Ca 2+ waves that create Ca 2+ gradients which play an important role in cell migration [1]. In addition, Ca 2+ waves also regulate the local transepithelial ion transport to maintain the spatial ion gradients across the epithelium [1]. We recently demonstrated in RPE that an easily induced and repeatable Ca 2+ wave could be produced by mechanical stimulation [5]. This provides an experimental way to study Ca 2+ activity in the epithelial monolayer.
In silico models of various cellular processes are becoming an increasingly important part of biological research, including drug discovery and toxicology studies. The importance of this was recently emphasized in a review of cardiotoxicity testing [10]. Computational models of Ca 2+ signaling, specifically, have been developed for many cell types including pancreatic and parotid acinar cells [11], astrocytes [12], and hepatocytes [13]. Epithelial Ca 2+ signaling, however, differs from other cell types because the epithelium forms a highly polarized cell monolayer that comprises organized apical and basal cell membranes. The epithelial cells are tightly connected with tight junctions and gap junctions between the cells [14]. At present, there are only a few epithelial Ca 2+ signaling models available, for example for the urothelial monolayer [15] and for the airway epithelium [16]. RPE has many unique functions compared to other epithelia as it supports the complex processes of vision. Indeed, in the treatment of many eye diseases, RPE is either the drug target or it hinders drug penetration and provides a barrier between most of the eye and the blood stream. Hence, computational models of the functions of RPE, including Ca 2+ dynamics, are well warranted.
The aim of this study, therefore, is to provide a deeper understanding of the study of Ca 2+ activity by introducing a detailed computational model of RPE Ca 2+ dynamics. The computational model described in this paper is based on our experimental data on a mechanically induced Ca 2+ wave in ARPE-19 cells, a commercial immortalized human RPE cell line that is widely used to assess RPE cell functions in vitro [17][18][19], regardless of its limitations in cellular morphology, organization and function [20].
The computational model is mostly based on the experimental data of Abu Khamidakh et al. 2013 [5]. In addition, the model comprises our new unpublished α-glycyrrhetinic acid (GA)-suramin-treated data. We constructed the model by combining previously published cell Ca 2+ dynamics model components of P 2 Y 2 receptors [21], inositol 1,4,5-trisphosphate (IP 3 ) receptors [22], ryanodine receptors [23], Ca 2+ pumps and gap junctions to a new model component of mechanical stretch. Furthermore, we connected the epithelial cells to each other in the model to enable the diffusion of the molecules and propagation of the stretch. We developed the model based on two experimental data sets: the GA-treated data, where gap junctions (GJs) were blocked by α-glycyrrhetinic acid and untreated control data, where GJs define the connections between the cells. The varying conditions the cells are exposed to due to the mechanical stimulation were modeled by defining three location-specific variables: stretch, extracellular ligand concentration, and IP 3 receptor phosphorylation rate. In addition, we validated the model by simulating the combined blocking effect of GJs and P 2 receptors by GA and suramin. This way, we obtained the first RPE Ca 2+ signaling model, and we could reveal a deeper understanding of Ca 2+ activity.

Experimental data
In this study, the experimental data of Abu Khamidakh et al. 2013 [5] 2+ imaging, by loading them with the Ca 2+ -sensitive dye fura-2-acetoxymethyl ester. Single cell mechanical stimulation, membrane perforation of one cell, was induced with a glass micropipette. The intracellular Ca 2+ concentration transient travelled over the ARPE-19 monolayer starting from the mechanically stimulated (MS) cell, and spreading to the neighboring (NB) cells (Fig 1). The NB cells immediately surrounding the MS cell were defined as the first NB cell layer (NB layer 1 = NB 1 ); cells immediately surrounding the first layer were defined as the second NB layer (NB layer 2 = NB 2 ) and so on. The ratio of the emitted fluorescence intensities resulting from excitation at 340 and 380 nm (F 340 /F 380 ) was determined for each cell after background correction. Normalized fluorescence (NF), which reflects the changes in intracellular Ca 2+ concentration, was then obtained by dividing the fluorescence value by the mean fluorescence value before the mechanical stimulation. [5] The experimental work produced data in NF units. The computational model, however, is presented in absolute calcium concentrations. Due to the lack of absolute reference we consider the model predictions only relative.
Three data sets were simulated with the model. Firstly, in the GA-treated data set the gap junctions (GJs) were blocked by α-glycyrrhetinic acid (GA) (Sigma-Aldrich, St. Louis, MO, USA). Secondly, the model was verified with an untreated control data set that was based on the previous model-only the GJ model component was added. Thirdly, the model was applied to predict a combined blocking effect of GA and P 2 receptor blocker suramin (Sigma-Aldrich) with GA-suramin-treated data set. Each data set was averaged from at least three separate experiments.
The experiments with GA-suramin-treated ARPE-19 cells were not included in the original paper of Abu Khamidakh et al. 2013 [5]. The experimental details concerning the ARPE-19 cells as well as the experimental solutions, infrastructure, and protocols are presented in [5] with the following exception: the cells were incubated in a solution containing 30μM GA (incubation time 30 min) and 50μM suramin (incubation time 25 min) prior to mechanical stimulation. To receive representative data for each NB layer, the raw data was averaged so that the NF graphs were aligned by the starting time of mechanical stimulation, and the mean values were calculated for each NB with a one-second sampling rate. This was previously done for the control data set [5], but the averaging was performed also here for the GA-treated and GAsuramin-treated data sets.
Indirect immunofluorescence staining. ARPE-19 cells (p. 24, 27, 44, three cover slips from each passage) were cultured on glass coverslips for two days. For immunofluorescence staining, the samples were washed three times with PBS and fixed for 15 min with 4% paraformaldehyde (pH 7.4; Sigma-Aldrich) at room temperature (RT). After three subsequent washes with PBS, the samples were permeabilizing by a 15 min incubation in 0.1% Triton X-100 in PBS (Sigma-Aldrich) at RT. This was followed again by three PBS washes, after which the samples were incubated with 3% bovine serum albumin (BSA; Sigma-Aldrich) at RT for 1 h. Primary antibody Zonula Occludens (ZO-1) 1:100 (33-9100, Life Technologies) was diluted in 3% BSA PBS and incubated for 1 h at RT. Samples were then washed four times with PBS, and Numbering of the cell layers. Schematic representation of the location of the mechanically stimulated (MS) cell with respect to the neighboring (NB) cell layers: NB 1 is the first layer which is in direct contact with the MS cell; NB 2 is the second layer which is in direct contact with NB 1 , and so forth. White line segment marks an apothem (a) of a hexagon. followed by 1h incubation at RT with secondary antibody donkey anti-mouse Alexa Fluor 568 (A10037, Life Technologies) diluted 1:400 in 3% BSA in PBS. The washes with PBS were repeated again and nuclei were stained with 4 0 , 6 0 diamidino-2-phenylidole (DAPI) included in the mounting medium (P36935, Life Technologies).
Confocal microscopy and image processing. Zeiss LSM780 LSCM on inverted Zeiss Cell Observer microscope (Zeiss, Jena, Germany) with Plan-Apochromat 63x/1.4 oil immersion objective was used for confocal microscopy. Voxel size was set to x = y = 66nm and z = 200nm, pixel stacks were set to 1024x1024, and approximately 50-80 slices were acquired with line average of 2. DAPI and Alexa-568 were excited with 405nm and 561nm lasers and detected with emission windows of 410-495nm and 570-642nm, respectively. The images saved in czi format were processed with ImageJ (Rasband, W.S., ImageJ, U. S. National Institutes of Health, Bethesda, Maryland, USA, http://imagej.nih.gov/ij/, 1997-2014.) and assembled using Adobe Photoshop CS6 (Adobe Systems, San Jose, USA).

Construction of the model
The Ca 2+ model was constructed by combining six subcellular model components that included the stretch component designed in this study and the P 2 Y 2 receptor models of Lemon et al. 2003[21], the IP 3 receptor type 3 (IP 3 R 3 ) of LeBeau et al. 1999 [22], and the ryanodine receptor (RyR) of Keizer & Levine 1996 [23]. The GJ model component connected the neighboring cells. These model components with corresponding numbering and their rationale, hypothesized Ca 2+ wave propagation mechanisms as well as model equations (see chapter Detailed model equations) that were used for the NB layers and data sets are summarized in Table 1. The basis for the mathematical implementation is presented in Fig 2 as a schematic model.

Parameters and parameterization
The model parameters are represented in Table 2 and the parameters specific for each NB layer in Table 3. Most of the parameters were adopted from the models of Lemon et al. 2003[21], LeBeau et al. 1999 [22], and Keizer & Levine 1996 [23]. Typically, the volumes of ARPE-19 cells [5,24,25] and RPE cells [26,27] are variable. The cell width was approximated to be 14μm from the corner-to-corner of a hexagon and the height was 12μm [5].The cytoplasmic volume was approximated to be about 70% of the total cell volume [28]. Thus, a cytoplasmic volume (v) of 1.07 10 -15 m 3 was used in the simulations. The initial values, the values at time of mechanical stimulation, were taken mostly from the model of Lemon et al. 2003 [21]. The initial value 0.12μM for intracellular Ca 2+ concentration ([Ca 2+ ] i ) is an arbitrary value approximating the baseline Ca 2+ concentration determined from GA-treated data set for NB 5 -NB 10 layers using Matlab SimBiology Toolbox.
The rest of the parameters were fitted with Matlab SimBiology Parameter Fit Task: First, the parameter values, excluding SSCC and GJ model components, were fitted with GA-treated data set in NB5 layer. This layer has in general the largest Ca 2+ response from those NB layers that do not experience any stretch due to mechanical stimulation, according to our assumption. Secondly, the SSCC model component parameters, excluding the location-specific stretch (θ) parameter (see below), were fitted with the same GA-treated data set in NB1 layer that is assumed to have the largest stretch. These values were then used in all simulations for all data sets and NB layers. For the control data set with gap junctions, all other parameters were kept unchanged but the GJ related diffusion parameters, D Ca 2þ , D IP 3 , In Ca 2þ and In IP 3 , were fitted using NB1 layer. As a boundary condition we assumed that there is no outflow of IP 3 and Ca 2+ outside the epithelium, thus Out IP 3 and Out Ca 2þ were assigned to be zero.

Location-dependent parameters
Three parameters were assumed to vary according to the location of the cell with respect to the MS cell: stretch (θ) activating the stretch-sensitive Ca 2+ channels (SSCCs), the extracellular ligand concentration ([L]) [6,16,29,30], and the phosphorylation rate of IP 3 R 3 (α 4 ) [22]. Ca 2+ concentration was modeled separately in each NB layer. The distance (x) defines the distance of the NB layer from the MS cell centre that was calculated using the idealized hexacon RPE cell architecture (Fig 1) as where a is an apothem of the hexagon, 6 is the number of corners in the hexagon, s = 7μm is the length of the hexagon side and n = 1, 2, 3. . .10 according to the NB layer numbering. The stretch component was present in cell layers NB1-NB4. Stretch (θ) was parameterized in the GA-treated data set separately for each NB layer. The obtained parameters resulted in an exponentially decaying function corresponding to the decay of an amplitude envelope of a damped wave in a membrane [31]. This function was then used for modeling the stretch where x is the distance from the MS cell (R 2 = 0.9878).  Ligand diffusion in the extracellular space is modelled according to thin film solution to Fick's diffusion law [32] as follows describing the ligand concentration (L) as a function of time (t) where L 0 is the initial bolus ligand concentration above the MS cell (at x = 0), D ATP is the diffusion coefficient for ATP, and x describes the NB layer distance from the central MS cell. IP 3 R 3 phosphorylation rate (α 4 ) used in Eq 21 was fitted separately for each NB layer in GA-treated and control data sets, which resulted in shallowly rising exponential functions with respect to the distance of the cell from the MS cell (x). The equation for GA-treated data set (R 2 = 0.9740) is and for control data set (R 2 = 0.9798) Similarly to θ and L, these functions were then used in simulations instead of values from separate fits.

Model simulations
The parameters were fitted with Matlab SimBiology (R2012a, The MathWorks, Natick, MA) to the experimental data using Parameter Fit task, where the maximum iterations was 100. The solver type was ode45 (Dormand-Prince) and the error model was constant error model. The time step in the simulations was set to Δt = 0.1 seconds.

Sensitivity analysis
Sensitivity analysis was performed to evaluate the uncertainty of selected parameters that were fitted in this study (parameters k IP 3 R 3 , k RyR , V Pump , K Pump , J Leak , In IP 3 , In Ca 2þ , D IP 3 and D IP 3 from Table 2) or behaved as location-specific parameters (parameters θ, L and α 4 from Table 3). Values of these parameters were changed -25%, -10%, 0%, +10% and +25% in the model including all the model components I-VI for the control data set. The influence of these parameter were studied for NB layers NB1, NB5 and NB10 concentrating on the following features of the Ca 2+ wave: peak amplitude, time to peak, Ca 2+ wave width at half maximum, and Ca 2+ concentration at the end of the Ca 2+ wave (at 90 seconds' time point).

Model prediction of drug effect: suramin
With the model, we investigated the mechanism by which suramin influences the Ca 2+ waves in ARPE-19 cells. First, we compared the peak amplitude, time to peak, Ca 2+ wave width at half maximum, and Ca 2+ concentration in the end of the Ca 2+ wave at 90 seconds' time point between two experimental data sets: GA-treated and GA-suramin-treated data sets. Second, we made sensitivity analysis about the behaviour of P 2 Y 2 receptor regulation parameters (K 1 , K 2 , k r , k p , k e , ξ), since suramin is a known unspecific antagonist of P 2 receptors. Suramin has also been suggested to disrupt the coupling between the receptor in the cell membrane and the Gprotein by blocking the association of the G-protein α and βγ subunits [33]. Hence, the G-protein cascade parameters k a , k d and δ were also evaluated. The sensitivity analysis was done in the model for GA-treated data set (including model components I-V) in NB1, NB5 and NB10 layers. The parameter values were changed in the model by -25% and +25% in order to compare the effects of parameter modifications to the observed differences in the experimental data between GA-treated and GA-suramin-treated data sets. All other parameters were kept unchanged. Third, based on the results of this approach, the model was fitted to the GA-suramintreated data set by refitting those P 2 Y 2 receptor and G-protein cascade parameters that were observed to change the Ca 2+ curve similarly to the differences seen in the experimental data between GA-treated and GA-suramin-treated data sets. This was done with Matlab SimBiology Parameter Fit task for each NB layer.

Detailed model equations
where the subscripts indicate the source of the flux: stretch-sensitive Ca 2+ channels (J SSCC ), inositol 1,4,5-trisphosphate (IP 3 ) receptor type 3 (J IP 3 R 3 ), and ryanodine receptor (J RyR ). J Pump combines the Ca 2+ pumping functions of sarco/endoplasmic reticulum ATPase (SERCA) and the plasma membrane Ca 2+ ATPase (PMCA). Leak Ca 2+ current (J Leak ) describes the total leakage from the extracellular space and the endoplasmic reticulum (ER) to the cytoplasm. J GJ;Ca 2þ is the Ca 2+ flux through gap junctions. I Stretch-sensitive Ca 2+ channels (SSCCs). Stretch-sensitive Ca 2+ channels (SSCCs) on the cell membrane are activated, when exposed to mechanical stimulation. Their closure is caused either by relaxation in the mechanical force or by their adaption to that mechanical force [34]. The SSCC model is described with Eqs 7-9. In this study, a model for SSCCs was developed according to the kinetic diagram shown in Fig 3, where C SSCC describes the proportion of the channels in the closed state. O SSCC is the proportion of SSCCs in the open state defined as where k f is the forward rate constant and k b is the backward rate constant. Ca 2+ flux via SSCCs (J SSCC ) is expressed as where k SSCC is the maximum Ca 2+ flux rate via SSCCs. Parameter θ is dimensionless, and describes the quantity of stretch induced at the time of mechanical stimulation, which then decreases with time according to a stretch-relaxation parameter k θ . II Purinergic receptor P 2 Y 2 . The agonist-induced activation of the second messenger system, here P 2 Y 2 , is represented by Eqs 10-16 [21].The kinetic diagram for the P 2 Y 2 receptor is presented in Fig 3. Some of the ligand-bound P 2 Y 2 receptors on the cell surface are phosphorylated irreversibly at rate k p , which causes desensitization of the receptors. Phosphorylated receptors are internalized at a rate k e , and these internalized receptors are then dephosphorylated and recycled back to the surface at rate k r . G-proteins can only be activated by the unphosphorylated P 2 Y 2 receptors [R S ] defined by where [R T ] denotes the total number of surface receptors, K 1 is the dissociation constant for unphosphorylated receptors, and [L] is the extracellular ligand concentration. The total number of phosphorylated surface receptors ½R S p is where K 2 is the dissociation constant for phosphorylated receptors. The binding of the ligand to the G-protein coupled receptor P 2 Y 2 results in a cascade of events leading to the activation of enzyme phospholipase C (PLC). This enzyme then hydrolyses the phosphatidylinositol 4,5-bisphosphate (PIP 2 ) to IP 3 . The activation rate (k a ) of the G-protein is proportional to two ratios: the ratio of the activities of the ligand unbound and bound receptor species (δ), and the ratio of the number of ligand bound receptors and the total number of receptors (p r ). Denoting the deactivation of G-protein to occur at a deactivation rate of k d , the equations for the amount of Gα Á GTP labeled as [G] as well as for the ratio p r can be expressed as and Equation for the concentration of IP 3 is where k deg is the degradation rate of IP 3 and J GJ;IP 3 is the IP 3 flux through gap junctions. The rate coefficient for PIP 2 hydrolysis (r h ) includes the effective signal gain parameter (α) and the dissociation constant for Ca 2+ binding to PLC (K 3 ) that can be expressed as Replenishment of PIP 2 is required for IP 3 production to be maintained over sustained periods of agonist stimulation. The equation for the number of PIP 2 molecules [PIP 2 ] is where r r represents the PIP 2 replenishment rate and [(PIP 2 ) T ] the total number of PIP 2 molecules. [21] III IP 3 receptor type 3 (IP 3 R 3 ). The IP 3 receptor type 3 (IP 3 R 3 ) function is represented by the Eqs 17-21 [22]. The kinetic diagram for IP 3 R 3 is shown in Fig 3. The IP 3 -induced release of Ca 2+ from the ER through IP 3 where k IP 3 R 3 is the maximum rate of Ca 2+ release, and IP 3 R 3

comprises four subunits that all must be in the open state (O) for the receptor to conduct. The steady-state proportion of open receptors (O) is
Where ϕ function controls the sensitivity of IP 3 R 3 to [IP 3 ], and it can be expressed as with rate coefficients k -1 , k 2 , k 3 , and k 5 being constants. Coefficient k 1 describes a rate for IP 3 R 3 transition from shut state (S) to open state (O) where constant α 1 is the maximum rate of S to O transition, and β 1 is the [Ca 2+ ] i at which the rate is half of its maximum. Coefficient k 4 expresses the rate for IP 3 R 3 from the first inactivated state (I 1 ) to the second inactivated state (I 2 ). It can be expressed as where the I 1 to I 2 transition is agonist specific and involves a phosphorylation of IP 3 R 3 by kinase activity. This is defined by parameter α 4 that denotes the maximum rate of I 1 to I 2 transition, while β 4 denotes the value of [IP 3 ] at which the rate is half maximal. [22] IV Ryanodine receptor (RyR). The ryanodine receptor (RyR) dynamics were modeled by Keizer & Levine 1996 [23] with Eqs 22-24. In Fig 3 the kinetic diagram for RyR is illustrated. The Ca 2+ release from the ER through RyR (J RyR ) is defined by the maximum RyR flux rate (k RyR ) multiplied by the open probability (P RyR ) as where and where w 1 is the RyR sensitivity function and K a , K b , and K c are dissociation constants. [23] V Sarco/endoplasmic reticulum Ca 2+ ATPase (SERCA) and plasma membrane Ca 2+ ATPase (PMCA). J Pump combines the pumping functions of sarco/endoplasmic reticulum Ca 2+ ATPase (SERCA) and plasma membrane Ca 2+ ATPase (PMCA) where V Pump indicates the maximum flux rate of the pumps and K Pump states the [Ca 2+ ] i for half-maximal pumping rate. VI Gap junctions (GJs). Gap junctions (GJs) and the Ca 2+ flux via GJs (J GJ;Ca 2þ ) are modeled as where n is the number of the NB layer. The NB layer n receives Ca 2+ from the previous NB layer n − 1 and delivers Ca 2+ to the next NB layer n + 1 according to the concentration gradient. Similarly, IP 3 flux through GJs (J GJ;IP 3 ) is modelled as where n is the number of the NB layer. D Ca 2þ is the diffusion coefficient for Ca 2+ and D IP 3 is the diffusion coefficient for IP 3 . These diffusion coefficients do not take into account the open probability, regulation, or density of the GJs as they describe the actual movement of Ca 2+ and IP 3 from one NB layer to the next NB layer. As an exception to other NB layers, the fluxes from MS cell to NB 1 layer are modelled by parameters In Ca 2þ and In IP 3 for Ca 2+ and IP 3 , respectively. Similarly, the fluxes from NB 10 layer to distant cell layers are modelled with parameters Out Ca 2þ and Out IP 3 . Parameter A describes the area of the cell membranes connecting the neighbouring NB layers in the monolayer. The value for A is received by multiplying the area of one hexagon side, that is the length of the hexagon side (l = 7μm) times the height of the cell (h = 12μm), by the number of hexagon edges between the two NB layers as where n is the number of the NB layer (n = 1, 2, 3. . ..10). Each NB layer has six cells with three connecting sides and (n-1) 6 cells with two connecting sides (see Fig 1). In other words, the area (A) increases with distance from the central MS cell.

Polarization of the ARPE-19 monolayer
Polarization of the ARPE-19 monolayer was demonstrated by immunolabeling the tight junctions in the monolayer. Confocal microscopy image (Fig 4) shows that within 2 days the ARPE-19 cells have formed a monolayer where ZO-1 is localised continuously in the junctions of the cells, forming a homogeneous network. This can be taken as an indication of the polarization of the epithelial cell culture [35].

Ca 2+ signal propagation mechanisms
The fittings of the model to the experimental data in the NB1-NB10 layers are illustrated in Fig  5A for the GA-treated data set and in Fig 5B for the control data set. The model simulations managed to catch very well the features of the experimental data in both data sets. In GA-treated data set (Fig 5A), the simulations closely followed the data in peak amplitude, time to peak, Ca 2+ wave width at half maximum and end Ca 2+ concentration in NB1-NB9 layers. In NB10 layer, however, time to peak was longer in the simulation results than in the data. In the control data set (Fig 5B), the Ca 2+ wave features differed slightly between the model and the data, but overall the curve shape of the model followed the data reasonably well. R 2 values describing the goodness of fit are presented in Table 4. In GA-treated data set and control data set R 2 values were higher than 0.8 in NB1-NB9 and lower than 0.8 in NB10. Hence, 90% of the fits in GAtreated data set and control data set resulted in R 2 > 0.8. The model includes the model components of SSCCs, P 2 Y 2 receptors, IP 3 R 3 s, RyRs, Ca 2+ pumps and GJs, and the parameters were either obtained from previous studies or defined in this study for ARPE-19. The basic fit was done in GA-treated data set for NB5, but the SSCC model component was fitted in NB1 (Table 2). Three location-specific parameters were defined in this study: stretch (θ), extracellular ligand concentration (L) and phosphorylation rate of IP 3 R 3 (α4) ( Table 3). The stretch (θ) and extracellular ligand concentration (L) decayed exponentially from NB1 towards the distant NB cell layers. The IP 3 R 3 phosphorylation rate regulated by the kinase activity (α 4 ) increased following a shallow exponential, almost linear function, from NB1 to NB10. The corresponding values of α 4 with the distance were lower in the control data set (Eq 5) than in the GA-treated data set (Eq 4) indicating a possible role of IP 3 receptor phosphorylation rate as a regulator of Ca 2+ signaling. The GJ model component was parameterized in control data set for NB1. GJs mediated the Ca 2+ signal by allowing the diffusion of Ca 2+ and IP 3 between adjacent cell layers so that the fluxes of these species decreased with distance from the MS cell due to the increasing area of the cell membranes connecting the NB layers ( Table 3).
The resulting model of mechanical stimulus induced Ca 2+ dynamics is: 1) Cells near the stimulus site conduct Ca 2+ through plasma membrane SSCCs, and gap junctions conduct the Ca 2+ and IP 3 between cells further away from stimulated cell. 2) The MS cell secretes one or several types of ligand to the extracellular space where the ligand diffusion mediates the Ca 2+ signal so that the ligand concentration decreases with distance.
3) The phosphorylation of the IP 3 receptor defines the cell's sensitivity to the extracellular ligand attenuating the Ca 2+ signal in the distance.

Results of the sensitivity analysis
The sensitivity of the four Ca 2+ wave features described in Materials and Methods was studied for a set of parameters that were fitted in this study for NB1, NB5 and NB10 layers (Fig 6). From the location-dependent parameters, θ, L and α 4 , Ca 2+ wave features were most sensitive to modifications in α 4 and the least sensitive to modifications in θ. Overall, decreasing the stretch parameter (θ) resulted in faster Ca 2+ waves in NB1 (Fig 6B). Increasing the extracellular   ligand concentration (L) in turn decreased the time to peak (Fig 6B) and increased the Ca 2+ wave width at half maximum ( Fig 6C). The effects of the changes in IP 3 receptor phosphorylation rate (α 4 ) on Ca 2+ wave features were complex: with decreasing α 4 Ca 2+ wave peak amplitude increased (Fig 6A), time to peak decreased (Fig 6B), the Ca 2+ wave width at half maximum increased or decreased depending on the NB layer (Fig 6C), and the end concentration increased (Fig 6D). The Ca 2+ wave features were insensitive to gap junction related parameters In Ca 2þ , D IP 3 and D Ca 2þ so that the tested modifications in their values resulted in less than 5% change in the features from the original conditions. Thus, they are not illustrated in Fig 6. However, increasing IP 3 input to NB1 via GJs (In IP 3 ) increased the Ca 2+ wave width at half maximum ( Fig 6C) and end concentration (Fig 6D), and the sensitivity was significantly higher in NB1 layer than in the more distant NB layers. In general, the sensitivity of the model to the changes in tested parameters depended on the NB layer and thus on the distance to the MS cell. Few parameters, however, were independent of the location (parameters k IP 3 R 3 , k RyR , V Pump , K Pump , J Leak ), but changes in their values affected significantly the investigated Ca 2+ wave features.

Possible suramin effect on attenuation of Ca 2+ waves
Similarly to the general sensitivity analysis, the four Ca 2+ wave features were compared between the experimental GA-treated and GA-suramin-treated data sets as well in NB1, NB5 and NB10 layers (Fig 7A). In the GA-suramin-treated data set, differences in peak amplitude were less than 10% compared to the GA-treated data set. However, time to peak decreased for NB1 and increased for NB5 and NB10 layers in the GA-suramin-treated data set. Furthermore, in this data set the Ca 2+ wave width at half maximum and the end Ca 2+ concentration were lower than in the GA-treated data set. Since suramin is a known P 2 receptor blocker, we studied the sensitivity of the model to P 2 Y 2 receptor parameters for the GA-treated data set (including model components I-V) by changing their values by ±25%. Similarly, G-protein cascade parameters were studied as well to consider the possible effect of suramin to disrupt the coupling between the receptor in the cell membrane and the G-protein. Our aim was to investigate the degree to which the observed differences in the Ca 2+ wave features between GA-treated and GA-suramin-treated data sets could be accounted for by the changes in these parameters. The sensitivity analysis revealed that the modifications in P 2 Y 2 unphosphorylated receptor dissociation constant (K 1 ) and P 2 Y 2 receptor phosphorylation rate (k p ) indeed induced changes that were similar to the experimental observations (see Fig 7A): increase in K 1 modified the time to peak (Fig 7B) and increase in k p narrowed the Ca 2+ wave width at half maximum ( Fig 7C). This would indicate disrupted ligand binding to the receptor or a higher phosphorylation rate of P 2 Y 2 receptors as well as a faster desensitization of the receptors after ligand binding. P 2 Y 2 receptor parameters K 2 , k r , k e , ξ, on the contrary, had a negligible influence on Ca 2+ wave behaviour: the modifications of these parameters by ±25% resulted only in less than 3% change on Ca 2+ wave features, as was the case also for G-protein cascade parameter δ. G-protein cascade parameters G-protein activation rate (k a ) ( Fig 7D) and G-protein deactivation rate (k d ) (Fig 7E) had diverse effects on the Ca 2+ wave features: for example decreasing k a and increasing k d narrowed the Ca 2+ wave width at half maximum in NB1 and NB5, but widened it in NB10. Thus, their behaviour did not follow the observations from the experimental data and therefore these factors were not considered to be responsible for the effects of suramin on the Ca 2+ wave. Fig 7F illustrates the fit of the model to the GA-suramin-treated data set after refitting the model parameters K 1 and k p . The values of these parameters ranged as follows: K 1 values decreased from 8.83μM in NB1 to 5.15μM in NB10 and k p values decreased from 0.19s -1 in NB1 to 0.05s -1 in NB10. Overall, the values of K 1 and k p were higher in the GA-suramin-treated data set than in the GA-treated data set. The simulated curves seem to fit well to the experimental data, except in NB1 layer at the end of the Ca 2+ wave. In GA-suramin-treated data set, R 2 values were higher than 0.8 in NB1-NB7 and lower than 0.8 in NB8-NB10 (Table 4). 70% of the fits in GA-suramin-treated data set resulted in R 2 > 0.8 indicating that the model explains only partially the combined effect of GA and suramin on Ca 2+ waves especially in the distant NB layers.

Discussion
The ARPE-19 cell line is an important biological model of human RPE despite its certain limitations [20]. This paper presents the first computational RPE model of Ca 2+ signaling using the experimental data measured from the ARPE-19 monolayer after mechanical stimulation. We aimed to create a model that combines the most important Ca 2+ signaling mechanisms in ARPE-19 cells so that the model can be used later in the development of more complicated RPE and epithelial models. Furthermore, the model was used to simulate and explain the Ca 2+ signaling of epithelia, especially RPE, taking into account the following factors: 1) cells are on the monolayer; 2) they are connected to each other by GJs permeating Ca 2+ and IP 3 , and 3) the cells are most probably experiencing different stretching and chemical conditions depending on their distance from the mechanical stimulation site. To the best of our knowledge, this is the first time as Ca 2+ signaling model has been implemented for the ARPE-19 monolayer. The model uses a set of location specific parameters including stretch, extracellular ligand concentration, and IP 3 R 3 phosphorylation rate as well as the Ca 2+ and IP 3 fluxes through GJs.

The identity of the extracellular ligand
The airway epithelium secretes the signal carriers ATP or UTP to the extracellular space in response to mechanical stimulation [16,36]. The connection of these ligands to Ca 2+ signaling as extracellular signal mediators has been mathematically modeled [16]. It is likely that a similar function can be linked to ARPE-19 or RPE, where the ligand interacts with the cell membrane P 2 Y 2 receptors. In our model, the ligand carried the signal in the extracellular space from the MS cell towards the distant NB cell layers after mechanical stimulation. According to our model, the extracellular ligand concentration decreased exponentially from NB1 towards NB10. We suggest, based on our modeling results, that the MS cell secretes ligand to the extracellular space. Epithelial cells such as ARPE-19 have been shown to secrete ATP under different stimuli [37,38]. On the other hand, the ligand degradation by ectonucleotidase activity [39] decrease the ligand concentration. The model predicts that the magnitude of the extracellular ligand concentration partly defines the nature of the cell response: higher and faster Ca 2+ waves were observed with higher ligand concentrations. The ligand concentration was derived from diffusion equation, and the obtained exponential decay function fitted well to the experimental data. Experimental studies show that the Ca 2+ wave peak amplitude value increases with increased ligand concentration in cultured human RPE [30] and in human airway epithelium [6]. Also, in the mathematical model of Warren et al. 2010 [16], it was observed that the time to peak for human airway epithelium decreased as the ligand concentration increased. These observations are in good agreement with our model.

The role of IP 3 receptor phosphorylation rate
The phosphorylation of the IP 3 receptor represents an important regulatory mechanism for Ca 2+ release [40][41][42]. It has been shown that the production of cyclic AMP (cAMP) through the activation of the adenylyl cyclase pathway leads to the activation of protein kinase A that phosphorylates IP 3 receptors [43].
Our simulation results show that the maximal phosphorylation rate of IP 3 R 3 (α 4 ) followed a shallow exponential, almost linear, increase from NB1 to NB10 in all three data sets. The parameter α 4 has previously been modeled as agonist specific only [22]. It is of note, however, that in addition to ATP or UTP and their interaction with P 2 Y 2 receptors, also other types of ligand-receptor interactions may occur. One plausible explanation could be that MS cell secretes different types of ligands, because its cell membrane was broken in mechanical stimulation. This would further lead to complex biological interactions at the cellular level, which is seen as a chance of this parameter with cell location. The need to model α 4 separately for the GA-treated data set and the control data set may be related to the functioning of the GJs, especially to their ability to alter ligand secretion in different cell types. Previous studies show that GJs participate in the regulation of the release of signaling molecules to the extracellular medium [44]. In astrocytes, as an example, GJs have been proposed to regulate the release of glutamate [45], an excitatory neurotransmitter and an important regulator of astrocyte Ca 2+ oscillations [46].
Overall, α 4 parameter may reflect a number of ligands and cell mechanisms not modelled in this nor other epithelial Ca 2+ models. The low α 4 values near the MS cell enable higher and faster Ca 2+ waves at corresponding ligand concentrations compared to the distal cell layers, where higher levels of kinase activity attenuate and slow down the signal. This aligns well with the literature. In RPE, the addition of 8-Br-cAMP counteracted the elevation of [Ca 2+ ] i induced by connective tissue growth factor (CTGF) [47], and the cell migration inhibitor adrenomedullin increased intracellular cAMP and decreased [Ca 2+ ] i [48]. The effect of the adenylyl cyclase pathway on IP 3 R kinetics has been ignored in most of the previously published Ca 2+ models e.g. [16,21,29], possibly because the kinase activity may not have been activated in those cell types or experimental conditions.

Gap junctions in Ca 2+ wave propagation
GJs connect the adjacent cells together and allow the diffusion of signaling molecules between them. The diffusion through GJs has previously been modeled, for example, in airway epithelium [16]. In our model, GJs carried the Ca 2+ signaling molecules between the NB layers based on the Ca 2+ and IP 3 concentration gradients, and permeated Ca 2+ and IP 3 selectively. As expected, NB layers near the MS cell were more sensitive to IP 3 input than the distant NB layers, and this was seen especially in the end Ca 2+ concentration at 90 seconds' time point.

Possible Ca 2+ wave attenuation mechanisms of suramin
In the GA-suramin-treated data set, the experimental data was reproduced in our model by increasing the unposphorylated receptor dissociation constant, which likely reflects disrupted ligand binding, and by increasing the phosphorylation rate of the P 2 Y 2 receptors to enhance their desensitization. This may indicate that suramin targets on P 2 Y 2 receptors as an unspecific P 2 receptor antagonist attenuating the Ca 2+ wave. This intriguing model hypothesis driven from the model results needs to be confirmed experimentally. It is worth noting, however, that suramin has also been considered to disrupt the coupling between the receptor in the cell membrane and the G-protein by blocking the association of the G-protein α and βγ subunits [33]. In our model, modifications in G-protein cascade parameters influenced the peak amplitude, time to peak and Ca 2+ wave width at half maximum. Despite the observed diversity in their effects between the NB layers, it is possible that suramin targets the G-protein cascade as well, by acting as an attenuator of the Ca 2+ wave.

Limitations of the model
Our work presents a computational model of epithelial Ca 2+ signaling based on experimental work on the ARPE-19 cell line. This cell line is used extensively as a model of RPE, although it differs from it to some extent. The limitations of ARPE-19 compared to native human RPE arise, for example, from cell organization and metabolism [20]. Importantly for our study, ARPE-19 cell line in our experimental setup lacked pigmentation which resulted in a lack of the large Ca 2+ stores, melanosomes, and needs to be taken into account when expanding our model to describe native RPE. In addition, we confirmed the polarity of the ARPE-19 monolayer with confocal microscopy. Trans-epithelial resistance (TER) that is a general measure of epithelial integrity was not measured due to technical challenges to perform the measurements on glass cover slips with our present equipment [49]. Nevertheless, the computational model created in this study describes the most important components of epithelial and ARPE-19 Ca 2+ activity. Thus it provides a good basis to address the native RPE in the future, even though it, being based on an in vitro model of RPE, needs to be considered only as a model. To improve the model further, experimental data and model implementations on certain additional Ca 2+ related mechanisms, such as P 2 X receptors [50], voltage-sensitive Ca 2+ channels [51] and Na + /Ca 2+ exchangers [52] would be well warranted. Finally, it is worth noting that the experimental work of Abu Khamidakh et al. 2013 [5] did not produce absolute Ca 2+ concentrations, and therefore our model also features only relative Ca 2+ activity.

Conclusions
A full mathematical understanding of RPE and epithelial Ca 2+ signaling would allow one to simulate cellular Ca 2+ responses under several physiological, pathological, and experimental conditions. Our present model represents significant progress towards this goal since it is able to reproduce the experimental data from an RPE type epithelium, ARPE-19 cell line, in different conditions, simulate several epithelial Ca 2+ signaling mechanisms, and predict drug responses in the epithelia. Our future work will include further development of the model especially focusing on the role of the voltage sensitive Ca 2+ channels in the RPE.