A depiction of the cell and corresponding calcium levels is given in Figure 1. The gating of the CRAC channels and to some extent the pump rate of the SERCA pumps depend on the calcium concentration inside the ER (Liu et al., 2010). The ER in turn is influenced by calcium entry through SERCA pumps, with which the cytosolic calcium concentration changes dynamically and which thus needs a separate description (Schmeitz et al., 2013). Taking the ER-PM junction into account, this results in three necessary calcium domains for a reliable model description. However, a real cell also contains additional important calcium pools such as the mitochondria or lysosomes, which are not considered at this level of abstraction. Additionally, the CRAC current in the SOCE model is assumed to be the only calcium influx into the cell, which however leaves out other calcium channels such as the TRP channels. The calcium pools are designed by the law of mass balance (Liu, 2012) which defines them by the sum of the calcium in- and outflows. All included calcium pools and flows as considered in this modeling approach are summarized in Figure 2.

Calcium enters the cytosol by diffusion from the junction, from the ER via the RyR and IP3R channels, and through an ER calcium leakage. The cytosol is cleared of calcium via the PMCA and the SERCA pumps, which are located in the plasma and ER membrane, respectively. The ER is filled with calcium via SERCA pumps, which are active in the junction and the cytosol, depending on the dynamics of local calcium concentrations. Therefore, two separate SERCA flows for the junction and the cytosol were implemented. The outflow from the ER is facilitated by IP3R channels, RyR channels and the ER leak. The majority of intracellular calcium is bound to calcium storage or buffer proteins (Liu et al., 2010; Islam, 2020). To account for buffered calcium, buffer terms were included for the cytosol and the ER, but not for the junction.

To verify whether the SOCE model actually describes and reproduces the calcium dynamics in A549 cancer cells, the model is compared with calcium measurements reported in the literature. Two studies were found that conducted suitable experiments on calcium channels important for this model (Hou et al., 2011; Shin et al., 2018). The experiment by Hou et al., 2011 was chosen as a reference since it directly shows the function of the CRAC channels and indirectly the function of IP3R channels, RyR channels and PMCA pumps. (see Figure 5C in Hou et al., 2011). The measurement depicts the relative fluorescence, which allows for conversion to the calcium levels (Padar et al., 2004). However, as the fluorescence parameters are not clearly defined in this work, a direct relationship and thus a quantitative comparison with the calcium levels is not possible. In detail, the measurement starts in a calcium-free environment, with a constant fluorescence intensity, indicating equilibrium of the cytosolic calcium level. Approximately 1 minute later, 2 µM Thapsigargin (TG) is added in order to block the SERCA pumps. This causes a net calcium efflux from the ER via IP3R and RyR channels as well as the ER leak, that raises the cytosolic calcium. The outflow continues until the ER is depleted within 2 minutes. At this point all STIM1 proteins should be bound to ORAI1, resulting in the opening of all available CRAC channels. In a next step, the cells are perfused with 2 mM Ca²⁺, which leads to a rapid and large rise of the calcium level as calcium flows into the cytosol, predominantely via CRAC channels. The steep rise is followed by a plateau phase, indicating the maximum calcium inflow in the cytosol through the fully activated CRAC channels. As the SERCA channels remain blocked, the ER is further depleted of calcium and the CRAC channels stay open. The calcium influx is then terminated by the addition of 100 µM 2-APB, a store-operated calcium channel blocker, causing the calcium levels to decline to the resting level, probably mediated by calcium extrusion pumps such as PMCA. To reproduce this experiment, the developed model was set to an arbitrary starting point with an external calcium concentration of zero. Depending on the parameter setting, the calcium concentration lasts some seconds to minutes until a steady state is reached. Addition of drugs and external calcium is assumed to occur at time points 66s for the TG, 231s for the perfusion with calcium and 429s for the CRAC channel blockage (Hou et al., 2011). The modulation of the corresponding channels depending on the effect of the drugs administered is then simulated by adjusting the model parameters accordingly.

Another suitable experiment is reported by Shin et al., 2018, where compound K was used to trigger ER Ca²⁺ release in A549 cells, showing the effect of the drug on cytosolic calcium levels, the results of which are depicted in Figure 4A in Shin et al., 2018. In the first experiment, compound K is added to the cell in order to trigger ER Ca²⁺ release, which is indicated by a quick and permanent increase in the cytosolic calcium concentration. In a second experiment the cells are incubated before addition of the channel blocker dantrolene and BAPTA-AM, an intracellular Ca²⁺ chelator, to block RyR and IP3R channels, which reduces the apoptotic effect and prevents increase of cytosolic Ca²⁺. While it is not possible to simulate the effect of compound K, it is possible to compute a forced calcium release through RyR channels by increasing the channels' flow rate in order to reproduce the measured data by model simulations (Shin et al., 2018).

Implementation of the SOCE sub-model and the combined A549-SOCE model using MATLAB/Simulink (R2019b, MathWorks Inc.) is described in the Supplementary Material Part C. As in the initial A549 model, the equations are solved using the forward Euler method (ODE1 in Simulink), with a fixed step size of 5・10⁻⁷ s. The combined model is run on a workstation (2–4 × 256 GB RAM, CPU 2.20 GHz (Intel(R) Xeon(R) Fold 5220R), 48 cores, 96 logical processors), with a runtime of approximately 3–4 h for a 10 s simulation. The final A549-SOCE model approach is shown in Figure 4.

Figures 5, 6 demonstrate the simulations for the settings of the SOCE concept (described in the Material and Methods section) based on the experiment by Hou et al., 2011 (cf. Figure 5C in Hou et al., 2011). Figure 5 shows a comparison of the cytosolic calcium concentration simulation versus the measured calcium fluorescence in the experiment. Note that the simulation yields an artificial steep calcium spike in the beginning, which reveals the transient response of the SOCE model setting to the calcium-free starting condition. Importantly, we emphasize that this transient response has no further impact on subsequent simulations. Once the initial phase is concluded, the simulation stabilizes, and the system reaches a physiologically relevant state. Despite all channels being active during this phase, the SERCA pumps are the most active, leading to a decrease of [Ca²⁺]_Cyt to a baseline of about 13 nM, while [Ca²⁺]_ER increases. The course of the buffered calcium inside the cytosol follows [Ca²⁺]_Cyt in shape, but with higher scaling. After 66 s, the flow rate of the SERCA pumps in the model is set to zero to simulate blockage by 2 μM TG application, resulting in a net outflow from the ER, and thus in a rapid decrease of [Ca²⁺]_ER and increase of [Ca²⁺]_Cyt. During the ER efflux the IP3R and RyR channels open and increase calcium release. After reaching a temporal maximum of 128 nM, [Ca²⁺]_Cyt declines as the PMCA responds to the rising calcium levels and clears the cell of the remaining calcium. In this phase, the CRAC channels are activated by the decreasing ER calcium and reach an open probability of 95%. The addition of extracellular calcium causes an immediate transient current through the open CRAC channels, which settles at 0.87 fA after a short settling phase. The influx is followed by PMCA outflow, resulting in an almost constant calcium concentration of 334 nM. The calcium concentration of the junctions follows the course of [Ca²⁺]_Cyt with roughly the same values up to this point. The CRAC channel influx causes [Ca²⁺]_Jun to reach 15 µM during the plateau phase. The endoplasmic calcium and the CRAC channel open probability remain unchanged. The subsequent blocking of the CRAC channels is achieved by setting the number of CRAC channels per junction to 0.001, causing [Ca²⁺]_Cyt and [Ca²⁺]_Jun to decrease rapidly. Figures 5D–F summarize the simulations of the relevant SOCE parameters [Ca²⁺]_ER, [Ca²⁺]_Jun, CRAC channel open probability, buffered Ca²⁺ inside the cytosol, and Ca²⁺ unbound calsequestrin.

Figure 6 shows the currents of the ion channels in the simulated SOCE experiment. During the initial resting phase, ER efflux provoked by IP3R and RyR currents is almost zero and calcium leaves the ER mainly through the implemented leakage current. During the SERCA inhibition after 66 s, the leakage decreases, but both, the IP3R and RyR currents, increase in a bell-shaped fashion. While these channels are active, the ER is completely emptied of calcium, as illustrated by the decrease of [Ca²⁺]_ER and increase in unbound calsequestrin in Figures 5B,E. The calcium-free calsequestrin is capped at 14 μM, which is assumed to be the total amount of protein and means that no more calcium is buffered in the ER. The PMCA current starts at about the same time as the ER depletion, but is larger and decays much slower, similar to [Ca²⁺]_Cyt. The addition of calcium leads to the opening of the CRAC and PMCA channels. Neglecting the spikes, i.e. the simulation artifacts caused by the direct supply of extracellular calcium, the CRAC and PMCA currents reach peak values of 0.9 pA. After CRAC blocking by 2-APB, only the PMCA current remains, which decreases rapidly. At the end of the simulation, all calcium compartments and buffers are almost empty and all calcium fluxes end.

To evaluate how different parameter settings affect the model’s behavior, various simulations were carried out (see Supplementary Materia Part D ). The simulations in Supplementary Figures SD1, A–C are shown for the cytosolic calcium concentration, as this parameter can be measured experimentally. Each sub plot in the figure yields a variation of a single parameter, while one graph (in blue) serves as the reference simulation with the basic settings used for Figure 5. These simulations support the basic functionality of the proposed model description, as the mechanism of the CRAC channels and the junction model work properly.

The simulation of the compound K experiment by Shin et al., 2018 as shown in Figure 8, starts in a steady state configuration. Hereto the SOCE model is run with an external calcium concentration of 1.8 mM for about 100 s, until steady state is reached. The resulting calcium concentrations with [Ca²⁺]_Cyt of 69 nM, [Ca²⁺]_ER at 660 μM, [Ca²⁺]_Jun at 0.9 µM and an open probability of the CRAC channels of 3% are shown in Figure 7. Note that the calcium concentrations depend on the parameter settings and thus result from the selected model parameterization. After reaching the baseline, the conductance of the RyR channel is increased by a factor of 10 at s 130, resulting in a steep increase of [Ca²⁺]_Cyt with a small overshoot that settles to a plateau after 20 s. The activation of RyR channels, rapidly depletes the ER of calcium and [Ca²⁺]_ER drops to 139 µM, as shown in Figure 8B. This, in turn, leads to an increased open probability of CRAC channels, which rises to 69% and lifts the junction calcium concentration [Ca²⁺]_Jun to 12 µM.

To compare simulations of the combined A549-SOCE model with the original A549 simulations, the combined model was run in steady state configuration for the cell cycle phases G0, G1, S and G2-M. Therefore, the membrane potential V_m was simulated by numerically solving Eq. 15 for V_m, starting at +5 mV until steady state was reached (t = 40 s, dt = 5.10^–7): C ∙ d V m d t = − I w h o l e _ c e l l (15)

Figures 9A–D show the calculated membrane potentials of the combined A549-SOCE model and of the original A549 in silico model for the cell cycle phases G0, G1, S and G2-M with the plasma membrane ion channel composition as denoted in the original A549 model version (Langthaler et al., 2021) for comparison, without taking into account potential alterations in ion channel expression within the ER or Ca²⁺ depletion from the ER during distinct phases of the cell cycle. Since the open probability of the CRAC channels is lower with the new Ca²⁺ formalism, the number of CRAC channels was set to 1500 instead of 200 in the original A549 model in order to achieve comparable CRAC currents. In addition, [Ca²⁺]_Cyt and [Ca²⁺]_Jun are also shown to examine the state of the SOCE model.

The membrane potentials as well as the response times of the combined model are roughly the same as those of the initial A549 model, as can be seen from the comparison of the membrane potential curves. For instance, the membrane potentials in G0 phase differ by less than 2 mV (V_m.original = −10.4 mV vs. V_m.calcium = −8.6 mV). A larger difference is seen in the S phase with V_S.diff = 5.9 mV (cf. Figure 9C). The combined model thus well reflects the oscillating characteristics of the membrane potential V_m during cell cycle progression (G0 phase: hyperpolarization, G1 phase: depolarization, S phase hyperpolarization, and G2-M phase depolarization)). Figure 10 outlines the rhythmic membrane potential change during cell cycle progression simulated with the initial and combined A549-SOCE model and relative changes of the calcium levels to the resting G0 phase. During the cell cycle phases, the cytosolic calcium levels fluctuate slightly between 0.182 and 0.204 µM and [Ca²⁺]_ER between 267.45 and 267.98 µM with the set parameters. The calcium levels of the junction are within 5.74 µM and 7.13 µM, which is higher than previously estimated by the initial A549 model. These high calcium levels result from increasing the maximum flow rate of the IP3R channels to 1.8 μM/s and the number of ER-PM junctions to 300, instead of the previously defined model parameters, to achieve currents in the expected sub nA range, as in the initial A549 model.

The SOCE model presented in this work allows the simulation and reproduction of the experimental data obtained by Hou et al., 2011 (cf. Figure 5C in Hou et al., 2011). However, the simulated calcium concentrations, in particular the initial cytosolic calcium concentration of [Ca²⁺]_Cyt = 13 nM could not yet be verified, as no experimental data in the absence of external calcium are available. In Padar et al. (2004) the baseline [Ca²⁺]_Cyt value of A549 cells is reported as 65 nM, presumably at an external calcium concentration of 1.2 mM. The [Ca²⁺]_Cyt levels reveal a relatively large difference but are in the same magnitude range, which indicates a reasonable model parameterization. The simulated ER calcium release matches the shape of SERCA blockage in the experiment in rise time and duration, but seems to have a steeper calcium decline. Since experimental values are missing, only a vague comparison of the simulated calcium plateau during the store-operated calcium entry is possible. Additionally, the ER Ca²⁺ release was compared with results of similar experiments of Padar et al., 2004, showing that the simulated data is consistent in terms of shape and time. However, using the stated relation and parameters for assessing the calcium concentration from the fluorescence response, the peak value of [Ca²⁺]_Cyt during ER depletion can be estimated to be about 360 nM, which is about three times the simulated peak value of 128 nM (Hou et al., 2011). This could indicate a higher calcium storage capacity of the ER in the real cell than assumed in the model simulations. The addition of external calcium leads to the activation of the store-operated calcium entry. The resulting [Ca²⁺]_Cyt during the SOCE plateau is 334 nM and thus very close to the expected [Ca²⁺]_Cyt of 300 nM determined by Hou et al., 2011. The other calcium compartments and the Ca²⁺ buffers equally demonstrate the expected behavior, but in contrast to [Ca²⁺]_Cyt cannot be further verified due to a lack of experimental data. The free calcium concentration of the ER is estimated to be in the range between 100 μM and 1 mM (Bong and Monteith, 2018; Rossi and Taylor, 2020), which matches the simulated values of almost 400 µM. However, the ER is completely emptied of free and bound Ca²⁺ during the ER depletion phase, which is reasonable from the simulation point of view, but may not be the case in a real cell considering the complex ER Ca²⁺ storage mechanisms (Islam, 2020). For [Ca²⁺]_Jun, estimates and simulations assume values of 1 μM–7 µM for the overall junction and up to 60 µM in the immediate vicinity of the plasma membrane (McIvor et al., 2018). Thus, the resulting [Ca²⁺]_Jun of 15 µM in the simulation is higher, but within the realm of possibility. By contrast, the values for the Ca²⁺ buffers cannot be substantiated as data on CSQ inside A549 cells is lacking and the cytosolic buffer is a simplified construct that only describes the combined response of several intracellular proteins and organelles. Overall, however, as could be shown, the model reliably simulates the effects of SERCA blockage by Thapsigargin and thus aligns well with known effects on calcium dynamics in cells from an electrophysiological point of view.

Steady state simulations, as shown in Figure 7, settle at [Ca²⁺]_Cyt of 69 nM under 1.8 mM extracellular calcium, which agrees well to the 64 nM at 1.2 mM extracellular calcium reported by Padar et al., 2004. With less extracellular Ca²⁺, the intracellular Ca²⁺ also becomes smaller, so it can be assumed that the model parameters used slightly underestimate [Ca²⁺]_Cyt. Nevertheless, the values are within the desired range, as [Ca²⁺]_ER settles at a high level of 594 μM, which indicates full Ca²⁺ stores.

The basic functionality of the model is underlined by additional parameter modulation shown in Supplementary Figures SD1, A–C ( Supplementary Material Part D ), as all parameter variations lead to the expected changes in simulated calcium concentrations. For example, increased IP3R and RyR conductivities in the simulation depict a faster Ca²⁺ efflux with a higher calcium peak, which is reasonable as these channels facilitate the ER efflux and more channels should result in a faster ion flux. Provided that no additional gating factors interfere, selective overexpression of these channels in a real cell should have a comparable effect which can therefore be simulated by the model approach. The number of junctions, ER leakage and diffusion distance constitute descriptions of cellular properties that are reasonable in a modeling context, but may not be directly applicable to a single feature in real cells due to abstraction. A valuable observation from these simulations is the relationship between [Ca²⁺]_Cyt and [Ca²⁺]_Jun, which can be seen by comparing Supplementary Figures SD1, A and Supplementary Figures SD1, B. The number of junctions and the PMCA rate variations have approximately the same effect on the junction and the cytosol. The diffusion distance, on the other hand, has a much greater effect on [Ca²⁺]_Jun than on [Ca²⁺]_Cyt and can therefore be used to adjust [Ca²⁺]_Jun independently of its environment. In addition, the ER simulation results are an important output of the model. Starting with the IP3R channels, it appears that they do not change the resting state of [Ca²⁺]_ER, but do affect the rate of ER efflux. Thus, the IP3R channel differs from RyR, which affects both the resting state [Ca²⁺]_ER and the rate of ER depletion. The leakage term has an unexpectedly large influence on the overall behavior, with too large leakage removing the Ca²⁺ from the ER and too small leakage significantly delaying the efflux. This shows that the ER leakage is an essential assumption for the model and is therefore probably also present in real cells. The increase of [Ca²⁺]_ER and the slower ER depletion with increasing amount of CSQ also represents a plausible simulation (see Supplementary Figures SD1, C) as CSQ represents a Ca²⁺ storage protein. The results of this variation suggests that the initial estimate of 14 µM CSQ was probably underestimated, as its influence on the model at concentrations below 140 µM appears to be small compared with other variables such as the cytosolic buffer capacity. A higher cytosolic buffer capacity, in turn, increases [Ca²⁺]_ER, which could be explained by the higher total Ca²⁺ available in the cell. Finally, since the SERCA pumps are blocked at the beginning of the experiment, the effect of SERCA variation can only be investigated in the initial resting phase. Here, an increased pump rate leads to an increase in [Ca²⁺]_ER, which also corresponds to the initial expectation.

Simulation of the experiment by Shin et al., closely matches the measured fluorescence response in time and shape, time and shape, thus reproducing the expected Ca²⁺ signaling chain for the CRAC channels. According to the experimental work, the simulation causes a permanent increase of [Ca²⁺]_Cyt due to the opening of RyR channels, as shown in Figure 8. Controlled Ca²⁺ depletion of the ER leads to the opening of CRAC channels and an increase of [Ca²⁺]_Jun, which subsequently influences the corresponding ion channels until a steady state is reached. While the calcium overshoot at the beginning, which is caused by the fast current through the RyR channels, could indicate a too pronounced opening rate in the model, all other simulation components behave as expected.

Considering that the presented model was not fitted to quality-assured experimental data, the simulated absolute calcium concentrations of the SOCE experiment are within physiologically reasonable ranges. Since the shape and scale of the simulation results can be widely varied via its parameters, as shown in the simulations, it is to be expected that the model can be fitted to suitable experimental data and thus a better agreement with the experiments can be achieved. The model obtains a reliable simulation of the [Ca²⁺] course in the junction and it also reproduces the forced ER Ca²⁺ release through RyR channels, providing a realistic estimate for the Ca²⁺ environment in the resting state of the cell. Overall, it can be concluded that at this stage the model can simulate the necessary calcium handling pathways that influence Ca²⁺ distribution of the ER-PM junction. The actual Ca²⁺ levels within the junction are consistent with estimates reported in the literature and can be easily adapted to the needs of the simulation. Therefore, the model can be used to make a reliable estimate of the Ca²⁺ concentration within the ER-PM junction.

The initial A549 model was designed assuming fixed values of [Ca²⁺]_Cyt = 64 nM and [Ca²⁺]_Jun = 4.68 µM. The cytosolic Ca²⁺ is important for the driving force of Ca²⁺ permeable channels, and Ca²⁺ within the junction is assumed to regulate the opening of KCa channels. In the combined A549-SOCE model these two Ca²⁺ parameters are dynamically calculated. Integration of the of the SOCE model into the A549 model required further adjustments in order to fulfil the initial assumptions, which are described below.

Simulation of the SOCE model resulted in a [Ca²⁺]_Cyt of 69 nM in the resting state and a much lower [Ca²⁺]_Jun of 0.9 µM. This makes sense from the Ca²⁺ perspective, as a low [Ca²⁺]_Cyt with non-zero external calcium suggests full Ca²⁺ stores and therefore low CRAC activation and [Ca²⁺]_Jun. However, the [Ca²⁺]_Jun proved to be too low to fully open the KCa channels, as assumed in the initial model, which requires [Ca²⁺]_Jun values above 5 µM. As a result, the current through the KCa channels was strongly reduced and the membrane potential depolarized considerably during first simulation trials. Hence, to ensure proper activation of KCa channels, the total calcium levels were increased by enhancing the maximum flow rate of IP3R channels. This has the dual effect of releasing more Ca²⁺ from the ER into the cytosol and increasing the influx of CRAC channels. As a second measure, the diffusion distance was adjusted to 350 nm for all simulations, increasing [Ca²⁺]_Jun independently of [Ca²⁺]_Cyt. The resulting calcium levels of [Ca²⁺]_Jun and [Ca²⁺]_Cyt allowed a better simulation of KCa channel gating, according to original model assumptions.

Figures 9A–D demonstrate the simulation of the membrane potential of the extendend model. By integrating the SOCE model, the membrane potentials V_m during the G0 phase differ by 2.9 mV and those during the G1 phase by 0.96 mV from the original A549 model. The ion currents of the displayed channels also have comparable values to the original simulations (data not shown), so that no significant change in channel functions can be observed. Only the CRAC current is slightly below the original range, which results from the lower opening probability of the CRAC channels with the new Ca²⁺ formalism, despite the increased channel number from 200 to 1500. The Ca²⁺ levels both in the junction and in the cytosol are significantly higher in both cell cycle phases, but are within the normal range. The Ca²⁺ transients take about 20 s to reach a steady state, which is slower than the voltage dynamics and expected compared with the general Ca²⁺ and SOCE dynamics in cancer cells.

Similar observations to the results of the G0 and G1 phases can also be made for the S phase and G2/M phase, which is shown in Figures 9C,D. The deviation in V_m for S phase of 5.9 mV is larger and thus slightly greater than expected, but could certainly be improved by using advanced computational optimization methods.

In summary, the simulation confirms that the SOCE model is compatible with the initial A549 ion current model and leads to an appropriate response of the calcium dependent ion channels to simulated Ca²⁺ levels. In contrast to the variation of V_m, the increasing cytosolic calcium concentration and the dynamic adaptation of the Ca²⁺ Nernst potential for the cytosol does not lead to a significant change in the driving force.

The new description of the CRAC channel provides a more profound model of CRAC gating, but its direct contribution to the whole-cell current is still small compared with the other Ca²⁺ channels. The important role of the CRAC channels for the model is regulation of the KCa channels. The initial model was designed with a high junctional calcium concentration, and thus, with almost fully open KCa3.1 channels, which was necessary to fit the experimental data. Although the high Ca²⁺ levels required for the simulation are within an acceptable range for cells, permanent activation of the CRAC channels is required to achieve the necessary calcium concentration for the KCa channels. However, the SOCE mechanism is only a temporary process limited to the time required to replenish the ER with Ca²⁺, and typical [Ca²⁺] levels in cells are much lower (Padar et al., 2004). In the case of lowering KCa3.1 activity, for example by reducing [Ca²⁺]_Jun due to the reduction of the diffusion barrierer, would lead to a depolarization of V_m above +10 mV. The problem could be explained by an insufficient KCa3.1 model that demonstrates a too low calcium dependence for the cell type. Another reason could be the data fitting of the initial model, as all cell cycle phases were fitted with the same [Ca²⁺] concentration. Similarly, the original assumption that the KCa channels are in close proximity to the CRAC channels may be incorrect. The KCa channels could instead be located near other calcium entry channels such as the TRP channels, which have a much higher conductance and could therefore generate higher local concentrations in their vicinity independent of the time limited SOCE mechanism. Nevertheless, more detailed data on the Ca²⁺ homeostasis are needed to allow adequate modeling of the actual KCa gating.

The entire SOCE concept of the developed model is based on several other model approaches. The regulatory mechanisms of the CRAC, SERCA, PMCA and IP3R channels are well known and have already been implemented with varying degrees of abstraction (Kowalewski et al., 2006; Means et al., 2006; Liu, 2012; Schmeitz et al., 2013; Eichinger et al., 2018). The novelty of this model lies in the unique combination of existing concepts and a higher level of detail. The most significant feature is the introduction of the ER-PM junction as a separate calcium domain, which has not been reported before in a similar form. Another sophisticated part is the addition of a kinetic state model for both the IP3R and RyR channels. While the IP3R channel has already been implemented as a kinetic state model (Means et al., 2006), this is typically an exception, and the RyR channel is generally neglected in SOCE models. Another notable sub-step is the Ca²⁺ and voltage dependency of the CRAC current. Although Ca²⁺ and voltage dependent CRAC current formalisms have been reported (Schmeitz et al., 2013; Eichinger et al., 2018), to the best of our knowledge, no model with the same level of detail was found in the literature reviewed.

An important geometric design decision is the diffusion barrier. It is generally assumed that the junction is directly surrounded by the cytosol on its side (Hogan, 2015; McIvor et al., 2018). However, the ER-PM junction is located in a microdomain with over 100 junctions, which means that the junctions that are not on the edge of the microdomain are not in direct contact with the cytosol. This would greatly enlarge the distance for the Ca²⁺ diffusion from the CRAC channels to the cytosol. If it is assumed that there is at least one junction with a diameter of 100 nm between the observed junction and the cytosol, the junctions are a few nm apart and the steady [Ca²⁺] decrease within the junction towards the edges is taken into account, then the diffusion distance for the junction would be around 120 nm–180 nm. This would mean that junctions in the center of the microdomain could have much larger diffusion distances and are almost isolated from the cytosol. With this assumption, a variable diffusion distance can be justified in the model, which can be set to determine [Ca²⁺]_Jun. The diffusion area, i.e., the lateral surface of the cylindrical junction, also influences the diffusion, but it depends on the size of the junction and therefore cannot be varied as easily. The other variables of the basic diffusion relationship are the dynamic calcium domains and the fixed diffusion constant, making the diffusion distance the parameter of choice for adjusting the ion flux between the junction and the cytosol.

The main Ca²⁺ influx into the junction occurs via the CRAC channels. The number of CRAC channels per junction could be varied to modulate [Ca²⁺]_Jun, but five channels is a number reported in the literature (McIvor et al., 2018) and was therefore retained for modeling. The gating formalism of the CRAC channels under the assumed simplification that the total amount of STIM1 and ORAI1 is simply 1 (µM) works as intended. The electrical driving force of the CRAC channels is regulated by the Ca²⁺ Nernst potential between the extracellular space and the junction. Simulation experiments during the development yielded much better results for a junction-based Nernst potential than for the [Ca²⁺]_Cyt -based one, as the rapid Ca²⁺ increase at the junction reduces the driving force and regulates the ion influx.

The chosen modeling approach for the SERCA is rather simple compared to other junction components, but it is sufficient to balance the calcium concentrations of the junction and cytosol and to fill the ER. No advanced kinetic state model has been found for the SERCA pump, and even the more detailed alternatives focus more on additional aspects such as the leakage (Means et al., 2006) or ATP dependency (Diederichs, 2008) rather than on the actual gating itself. The transport or reaction rate of the published SERCA models is typically only controlled by a single dissociation constant. The dissociation constant in this model does not represent real protein affinities, but was determined empirically and set to a value that yielded good results. Deviations from the real SERCA can be explained by the fact that there are different isoforms of SERCA pumps present in the cells (EMBL-EBI, 2023) and the formalism used represents an approximation of the combined behavior of different pumps with different Ca²⁺ affinities. Almost exactly the same observations and restrictions are made for the PMCA pump, with the difference that it only affects the cytosol. These pumps have the same basic functionality, are modeled by the same basic concept and give similarly satisfying results.

In contrast to SERCA and PMCA pumps, IP3R channels have already been described by various kinetic state models (De Young and Keizer, 1992; Sneyd and Dufour, 2002; Means et al., 2006). The modeling approach used provides expected results that fit well with the supposed SOCE mechanism. One limitation in modeling this channel is its dependence on the second messenger IP3. The highly simplified mathematical description of the IP3 dynamics used in this work does not reflect the biological processes that lead, for example, to an instant activation of the IP3R channel, but does allow a basic simulation of IP3R activation. Although such functionalities are beyond the scope of this first modeling approach, it must be considered a limitation. Similar model limitations also apply to the RyR channels, which are usually neglected in computational approaches. Although RyR channels play a minor role in ER Ca²⁺ release compared with IP3R, they are functionally expressed in A549 cells and can cause a significant ER Ca²⁺ efflux independently of IP3R (Shin et al., 2018). The addition of RyR channels therefore expands the possible simulation scenarios and represents the SOCE mechanism in more detail. These advantages justify the comparatively minor disadvantage of the additional computational effort in this stage.

The gating of the RyR model requires a description of the Ca²⁺ buffer CSQ, which leads to a further structural limitation of the model due to the introduction of buffer equations. The formalism for CSQ is appropriate and sufficient for RyR channel gating, but does not represent the buffered Ca²⁺ in the ER. Based on the results discussed it can be assumed that the initial amount of CSQ is probably too low and thus insufficient Ca²⁺ binds. It is also important to note that CSQ is not the major Ca²⁺ binding protein in the ER, and probably attributed to calreticulin (Means et al., 2006; Schmeitz et al., 2013), as a result of which the model does not simulate the buffered Ca²⁺ in the ER. In addition, the dimensioning of the entire ER is based on an assumption and the actual size of the ER could differ. Although the results for the Ca²⁺ inside the ER appear to be in the correct range and follow reasonable trends, these have to be interpreted with caution as they may deviate strongly from the actual Ca²⁺ environment of the real cell. The description of the buffer dynamics of the cytosol, on the other hand, represents a subsumption of the different Ca²⁺ binding proteins, which fulfills its purpose of shaping the Ca²⁺ dynamics to fit the reported experimental data but, like the CSQ description, does not accurately represent the amount of buffered Ca²⁺ of a real cell.

The TRP channels in the original A549 version are modeled as constant conductivities and therefore constantly open. The high conductivity causes large Ca²⁺ fluxes, which lead to a permanent overfilling of the SOCE model with Ca²⁺. For this reason, although the TRP channels significantly modulate the calcium concentration of the cells, they were excluded from the current SOCE model.

The main important limitation is the lack of experimental data for parameter fitting. Although the presented model can reproduce the SOCE mechanisms of the A549 cell line quite accurately in terms of time and shape of fluorescence intensity measurements, the calcium concentrations cannot be verified, which limits its predictive ability and reliability.

To overcome existing limitations and improve accuracy, parameter fitting through computational optimization and a reparameterization of the A549-SOCE model is an essential next step, provided that data are available that allow quantitative simulations of the SOCE mechanism. For instance, recent studies have yielded valuable insights and successfully measured the free calcium concentration inside the ER lumen (Rossi and Taylor, 2020), which can serve as a reference for refining the ER description. In general, incorporating additional data concerning the Ca²⁺ homeostasis of A549 cells, such as resting [Ca²⁺] levels at different stages of the cell cycle, would certainly improve the model’s performance. Moreover, there’s a need to gain deeper knowledge of cell specific ion channel dynamics, such as TRP6, TRPV3, CRAC, and KCa channels. Developing new hidden Markov models or kinetic state models for these channels could significantly refine the SOCE model, potentially improving its accuracy in predicting parameters such as membrane potential (V_m) and [Ca²⁺]_Jun.

The presented SOCE model simulates the dynamics of the CRAC channel and the Ca²⁺ concentration specifically within the ER-PM junction. However, despite considering other Ca²⁺ domains, it is not intended as a general exhaustive representation of calcium dynamics encompassing all facets within the cell. Developing a comprehensive calcium model would necessitate a more profound description of the major Ca²⁺ buffer proteins in the ER and cytosol as implemented, e.g. in Means et al., 2006. The addition of other Ca²⁺ permeable channels or pumps, such as the Na⁺/Ca²⁺ exchanger would equally improve the accuracy of the model. A more sophisticated model should also include the mitochondria, which are the second largest intracellular Ca²⁺ store and, particularly important for the Ca²⁺ homeostasis. In addition, a detailed description of the IP3 cycle for IP3 mediated Ca²⁺ release in the ER or the possibility of simulating Ca²⁺ oscillations could be implemented. Some of these examples have already been described and implemented elsewhere and should be compatible with the developed SOCE model (Marhl et al., 2008; Liu, 2012; Schmeitz et al., 2013).

The presented model focusing only on electrophysiological processes, sets the stage for future explorations. The path toward enhancing the A549-SOCE model is multi-layered and involves optimizing the parameters, incorporating experimental data and validation, and, ultimately, extending the model to the downstream signaling pathways. This iterative process of refinement is fundamental in bridging the gap between theoretical modeling and biological reality, thereby laying the basis for novel, comprehensive tools in cancer research.