1 Introduction

The endocrine pancreas comprises about 1 million islets scattered throughout the pancreas, covering up to 1–2% area. Each islet contains \(\sim\)1000 \(\beta\)-cells surrounded by the \(\alpha\), \(\delta\) and PP or F-cells responsible for the secretion of glucagon, somatostatin and pancreatic polypeptide, respectively Ferrannini (2010). These cells are excitable and use extracellular spaces and gap junctions to communicate with each other (Escher and Long 2013; Zimny and Blackard 1975). The pancreatic \(\beta\)-cells are the only source that secretes insulin in the human body responding to elevation in plasma glucose concentration. The dysfunction in these cells is one of the major causes of diabetes (Rorsman 2005; Jaberi-Douraki et al. 2015; Portuesi et al. 2013). Insulin has a significant role in metabolizing carbohydrates and fats, maintaining glucose homeostasis and proliferation, etc. Insulin is a small protein made up of two polypeptide chains. The threshold glucose concentration for insulin secretion is 3.3 mmol/l (Aronoff et al. 2004). The glucose stimulates insulin secretion, which is then involved in the metabolic pathway of sugar and elevates the intracellular ATP and ADP ratio. Due to this, the ATP sensitive potassium channel (\(K_{ATP}\)) present in the plasma membrane gets closed and results in depolarisation of the cell, which in turn causes an elevation in calcium influx by the voltage-dependent \(Ca^{2+}\) channels (VGCC), thereby increasing the cytoplasmic calcium. The cytoplasmic calcium acts as a second messenger in \(\beta\)-cells and regulates processes like ATP production, transcription, metabolism, insulin, amylin and C-peptide secretion (Idevall-Hagren and Tengholm 2020). The increase in cytosolic calcium activates insulin exocytosis in \(\beta\)-cells Pedersen (2009). Also, the endoplasmic reticulum (ER) is involved in regulating glucose levels in \(\beta\)-cells as ER acts as a dynamic storage facility for calcium (Idevall-Hagren and Tengholm 2020). The calcium release from ER is primarily caused by \(IP_{3}\)-receptor mechanisms. \(IP_{3}\) activates the \(IP_{3}\)-receptor, indicated by \(IP_{3}R\). Due to the leak attached to the ER membrane, free \(Ca^{2+}\) ions leak into the cytosol from ER, which breaks phosphatidylinositol-4,5 bisphosphate (\(PIP_{2}\)) to release \(IP_{3}\), which then binds with \(IP_{3}\) receptor (\(IP_{3}R\)) and opens the receptor gate to facilitate the diffusion of free \(Ca^{2+}\) ions from ER to the cytosol.

The \(Ca^{2+}\)-ATPases SERCA pump mediates the uptake of \(Ca^{2+}\) in \(\beta\)-cells to maintain the high [\(Ca^{2+}\)]\(_{ER}\), which is necessary for the survival and functions of \(\beta\)-cells. Although the \(Ca^{2+}\) influx and removal are maintained by channels and pumps inside the cell, there are other factors that help in maintaining the \(Ca^{2+}\) level inside the cell, known as buffers. As a high level of free calcium concentration is toxic for the cell, the buffers work as the binding proteins for calcium (Dupont et al. 2016). EGTA, BAPTA, Sorcin and Persinilin-1 are the major buffers found in \(\beta\)-cells (Idevall-Hagren and Tengholm (2020)).

A number of studies by theoretical and experimental approaches have been noticed working hand in hand to deduce new results and contribute to the progress of research on various cells. However, such types of studies have not been reported for \(\beta\)-cells. Various analytical and numerical models of calcium dynamics in different cells like myocytes (Pathak and Adlakha 2015a, b), astrocytes (Jha et al. 2012, 2013, 2020), oocytes (Naik and Pardasani 2015, 2013), hepatocytes (Jagtap and Adlakha 2018), acinar cells (Manhas and Pardasani 2014; Manhas et al. 2014), neurons Tewari and Pardasani (2009), fibroblast cells (Kotwani et al. 2012), etc. are reported in the literature.

Tewari and Pardasani (2009) analyzed the effect of sodium influx on the concentration of calcium in neurons. Manhas and Pardasani (2014) introduced a simple model for calcium changes in the pancreatic acinar cells involving the \(IP_{3}R\) and RyR (ryanodine receptor). Manhas et al. (2014) also performed numerical simulations to study complex calcium oscillations in acinar cells involving calcium modulation of \(IP_{3}\) levels by feedback regulatory mechanisms of degradation and production. Pathak and Adlakha (2015a, b) have constructed a mathematical model to analyze the calcium concentration levels in myocytes by incorporating calcium diffusion, excess buffers, pump and leak for a one-dimensional spatio-temporal case. Using the finite volume method, Jha et al. (2012) constructed a mathematical model to explore the relationships among the parameters impacting \(Ca^{2+}\) dynamics, such as source influx and diffusion coefficient in astrocytes. They also studied the effect of voltage-dependent \(Ca^{2+}\) channels (VGCC) on their model for two and three-dimensional space and obtained results by using the finite element method Jha et al. (2013, 2020). Kotwani et al. (2012) proposed a model to examine calcium ion concentration changes in a fibroblast cell involving excess buffers.Naik and Pardasani (2015, 2013) developed a one-dimensional mathematical model of spatio-temporal calcium variations in oocytes by incorporating the RyR, VGCC and diffusion coefficients. Further, they extended their model to higher dimensions to study the calcium regulation in oocyte cells (Naik and Pardasani 2017). Recently, the calcium distribution in T-lymphocytes has also been studied with the help of different dimensional models by incorporating the various parameters (Naik 2020; Kumar et al. 2018). Jagtap and Adlakha (2018) studied the \(Ca^{2+}\) variations in the presence of buffers and fluxes for two dimensions with the help of the finite volume method in a hepatocyte cell.

Wagner and Keizer (1994) analyzed the effect of buffers on the movement of \(Ca^{2+}\) in neurons from internal stores, and they also reported the effects of \(IP_{3}\), rapid stationary, and mobile \(Ca^{2+}\) buffers. Jafri and Keizer (1995) explored the impact of the mobile and stationary buffers, diffusion and the endoplasmic reticulum on \(IP_{3}\)-induced \(Ca^{2+}\) wave propagation. Various models based on the experimental data obtained from the rodent are reported for electrical activity-dependent \(Ca^{2+}\) signaling in a \(\beta\)-cell (Pedersen 2009; Fridlyand et al. 2003, 2009; Varadi et al. 1995; Sabatini et al. 2019; Keizer and Magnus 1989; Chay 1986). Chay (1997) elucidated the importance of the extracellular \(Ca^{2+}\) concentration regulating the luminal and intracellular \(Ca^{2+}\) concentrations in the intracellular \(Ca^{2+}\) stores and electrical signaling in a \(\beta\)-cell. Idevall-Hagren and Tengholm (2020) discussed the effect of glucose on cytosolic \(Ca^{2+}\) in a \(\beta\)-cell for generating signaling patterns in spatial and temporal directions involved in insulin secretion. They also discussed the coordination mechanism of \(Ca^{2+}\) signals among \(\beta\)-cells to understand the dysfunctions in Type-2 diabetes.

Keizer and Magnus (1989) constructed a model for the electro-physiological mechanisms of \(\beta\)-cells involving voltage-dependent outward potassium and inward calcium current mechanisms.

Fridlyand et al. (2007) developed mathematical models for cAMP dynamics regulation and fast calcium dynamics depending upon electro-physiological events for pancreatic \(\beta\)-cells. They also coupled mathematical equations for the critical ionic currents and calcium homeostasis in a \(\beta\)-cell and showed that changes in \(Ca^{2+}\) concentration can modulate inositol lipid-specific phospholipase activity, which in turn can stimulate \(IP_{3}\)-production (Fridlyand et al. 2013).Gromada et al. (1996) performed an experiment on insulin-secreting \(\beta\)TC3 cells and conclude that [\(Ca^{2+}\)]\(_i\) has an important role in controlling \(IP_{3}\)-production in \(\beta\)TC3 cells.

Seino et al. (2011) presented a theoretical study on how impaired insulin secretion can cause obesity and diabetes.

Tamarina et al. (2005) performed an experiment to study the relationship between \(IP_{3}\) and calcium oscillations in \(\beta\)-cells and concluded that \(IP_{3}\) oscillates in parallel with glucose-stimulated intracellular calcium.

Pedersen et al. (2005) formulated a mathematical model for insulin secretion produced due to calcium-driven electrical oscillations in the combination of glycolysis oscillations.

Rorsman et al. (2012) presented a model for reciprocal regulation mechanisms for glucagon and insulin secretion stimulated by plasma glucose for \(\alpha\) and \(\beta\)-cells. They also concluded that variation in calcium signaling can disturb the functions and performance of \(\beta\)-cells, leading to the development of Type-2 diabetes.

Few experimental studies have also been conducted on human islets and claim the significant role of calcium in the human body. Braun et al. (2008) identified the VGCC expressed in human \(\beta\)-cells obtained from non-diabetic donors and characterized their role in glucose-induced insulin secretion. They concluded that the \(Ca^{2+}\) primarily triggers the exocytosis of insulin-containing granules. Misler et al. (1992) tested an electrophysiologically-based hypothesis of the coupling of glucose metabolism to insulin secretion in human \(\beta\)-cells.

The review of the literature gives us an idea that most of the reported models have been developed for electrical activity-based calcium dynamics in \(\beta\)-cells. As far as, we are aware that many attempts have been made to study the intracellular calcium regulatory mechanism in other cells but not in \(\beta\)-cell. However only few attempts have been made to study the temporal dynamics of calcium in a \(\beta\)-cell involving SERCA pump and ER-leak. To the best of our knowledge, no model of spatio-temporal calcium dynamics regulating production and degradation of \(IP_{3}\), production of ATP and insulin secretion in a \(\beta\)-cell is reported in the literature. Thus there is a need to develop a model to investigate these impacts of calcium dynamics on other various processes like production and degradation of \(IP_{3}\), production of ATP and insulin secretion in a \(\beta\)-cell under various conditions of health and disease. The development of such models will aid in exploring dysregulatory impacts of calcium dynamics in a \(\beta\)-cell which can cause various metabolic disorders, dysfunction of multiple processes and their implication in various diseases like diabetes, obesity, etc. Therefore, the primary goal of the present study is to develop a spatio-temporal model for calcium dynamics which can further regulate the \(IP_{3}\), ATP and insulin production in a \(\beta\)-cell. Here, a mathematical model of calcium dynamics with one-way feedback of calcium in the regulation of \(IP_{3}\), ATP and insulin production has been developed. ATP production depends on both calcium and \(IP_{3}\) (Stamatakis and Mantzaris 2006). \(IP_{3}\) (Hirose et al. 1999; Manhas and Pardasani 2014; Politi et al. 2006) and insulin (Pedersen et al. 2005) production depends on calcium concentration in a \(\beta\)-cell. Thus, calcium also provides indirect feedback through \(IP_{3}\) in ATP-production. The \(\beta\)-cell plays an important role in glucose metabolism. Therefore, proper functioning of the \(\beta\)-cell is important for this glucose metabolism. The calcium dynamics has an important role in maintaining the structure and functions of the \(\beta\)-cell. Apart from this, \(IP_{3}\), ATP and insulin also have roles in the functioning of the \(\beta\)-cell. The focus of the proposed study is on the regulatory role of the calcium dynamics in \(IP_{3}\), ATP and insulin production in a \(\beta\)-cell.

2 Mathematical Model

The modeling investigations have begun by outlining the processes by which calcium is released into the cytosol. After establishing a quantitative understanding of the behavior of this calcium subsystem, the single-cell model is being extended in the present study to explore the \(IP_3\), ATP and insulin levels. The calcium dynamics in a \(\beta\)-cell that involves an excess buffer is expressed by the following partial differential equation (Smith et al. 1996; Crank 1979):

$$\begin{aligned} \frac{\partial [Ca^{2+}]}{\partial t}=D_{Ca}\nabla ^{2}[Ca^{2+}]-k_{i} ^{+}[B_{i}]_{\infty }([Ca^{2+}]-[Ca^{2+}]_{\infty })+J_{leak}-J_{SERCA}, \end{aligned}$$
(1)

where the first term on the right is the diffusion term and the second term is the reaction term. The symbols used in the above equation are defined below,

      \([Ca^{2+}]\): Cytosolic calcium concentration;

      \([Ca^{2+}]_{\infty }\): Free cytosolic \(Ca^{2+}\) concentration at rest;

      \([B_{i}]_{\infty }\): EGTA buffer concentration;

      \(D_{Ca}\): Diffusion coefficient for calcium;

      \(k_{i} ^{+}\): Rate of association for EGTA buffer.

$$\begin{aligned} J_{leak}=P_{ER}([Ca^{2+}]_{ER}-[Ca^{2+}]), \end{aligned}$$
(2)

where \(P_{ER}\) is ER-leak permeability and \([Ca^{2+}]_{ER}\) denotes the free calcium concentration in ER (Fridlyand et al. 2003).

$$\begin{aligned} J _{SERCA}=\frac{P^{max}_{SERCA}[Ca^{2+}]^2}{k_{pump} ^{2}+[Ca^{2+}]^2}, \end{aligned}$$
(3)

where \(P _{SERCA}^{max}\) is the maximum pumping rate of the SERCA and \(k_{pump}\) is the half-maximum pump activity of the SERCA pump (Dupont et al. 2016).

Production of \(IP_{3}\) depends on the activation of phospholipase-C (PLC), which also relies on the \(Ca^{2+}\) feedback control of cytoplasmic concentration (Hirose et al. 1999). \(IP_{3}\) degradation is also \(Ca^{2+}\)-dependent (Manhas and Pardasani 2014; Politi et al. 2006). Inositol(1,4,5) trisphosphate 3-kinase removes \(IP_{3}\) from the cytosol through phosphorylation. The temporal variation in \(IP_{3}\) concentration is defined as:

$$\begin{aligned} \frac{d [IP_3]}{dt}= {IP_{3}}_{production}-{IP_{3}}_{degradation}, \end{aligned}$$
(4)
$$\begin{aligned} {IP_{3}}_{production}=k_{prod}\frac{[Ca^{2+}]^2}{K_{prod}^2+[Ca^{2+}]^2}, \end{aligned}$$
(5)
$$\begin{aligned} {IP_{3}}_{degradation}=k_{deg}\frac{[Ca^{2+}]^2}{K_{deg}^2+[Ca^{2+}]^2}[IP_3]_{cyt}, \end{aligned}$$
(6)

where \(k_{prod}\): Maximal production rate of PLC isoforms; \(K_{prod}\): Sensitivity of PLC to calcium; \(k_{deg}\): Phosphorylation rate constant; \(K_{deg}\): Half saturation constant of \(IP_{3}\) kinase; \([IP_3]_{cyt}\): Cytosolic \(IP_3\) concentration at equilibrium (Ehrlich 2000).

ATP has a crucial role in stimulus secretion coupling in pancreatic \(\beta\)-cells (Ehrlich 2000; Idevall-Hagren and Tengholm 2020). In the literature, there have been two main hypotheses put forth. The first theory explains the \(Ca^{2+}\)-dependent ATP release described by Cotrina et al. (1998) and the \(IP_3\)-dependent ATP release described by Wang et al. (2000). The \(Ca^{2+}\)-dependent ATP production is given as (Stamatakis and Mantzaris 2006; Rorsman 2005):

$$\begin{aligned} ATP(Ca^{2+})=\omega ATP_{prod}^{max}\left( \frac{\frac{F_{0}}{F_{0}-1}-2(\frac{[Ca^{2+}]}{[Ca^{2+}]_{max}})}{\frac{1}{F_{0}-1} -(\frac{[Ca^{2+}]}{[Ca^{2+}]_{max}})^2}\right) . \end{aligned}$$
(7)

The \(IP_{3}\)-dependent ATP production is given as Stamatakis and Mantzaris (2006):

$$\begin{aligned} ATP(IP_{3})= \omega ATP_{prod}^{max}\left( \frac{\frac{G_{0}}{G_{0}-1}-2(\frac{[IP_{3}]}{[IP_{3}]_{max}})}{\frac{1}{G_{0}-1} -(\frac{[IP_{3}]}{[IP_{3}]_{max}})^2}\right) . \end{aligned}$$
(8)

The above \(Ca^{2+}\) and \(IP_{3}\)-dependent ATP production functions are constructed in such a way that \(ATP(Ca^{2+})\) (\(ATP(IP_{3})\)) reaches a maximum when [\(Ca^{2+}\)]= [\(Ca^{2+}\)]\(_{max}\) (\([IP_{3}]\)=\([IP_{3}]_{max}\)) and is equal to \(F_0\) (\(G_0\)) for [\(Ca^{2+}\)]=0 \(\mu M\) (\([IP_{3}]\)=0 \(\mu M\)). A small, non-zero value has been considered for \(F_0\) and \(G_0\) to describe basal ATP release which is 0.05.

As shown in the literature review (Pedersen et al. 2005; Rorsman et al. 2012), insulin secretion depends on calcium concentration. Finally, to describe the insulin secretion rate, the model has been considered as given by Pedersen et al. (2005):

$$\begin{aligned} \frac{d[I]}{dt}=\frac{I_{\infty }(Ca)-[I]}{\tau _{a}}, \end{aligned}$$
(9)

where \(I_{\infty }(Ca)\) denotes equilibrium secretion rate and is defined as:

$$\begin{aligned} I_{\infty }(Ca)= \left\{ \begin{array}{ll} I_{slope}([Ca^{2+}]-Ca_{null})\quad for &{} \quad [Ca^{2+}] \ge Ca_{null} \\ 0 \quad \quad for &{} \quad [Ca^{2+}] < Ca_{null} \end{array}, \right. \end{aligned}$$

where \(I_{slope}\) measures the \(Ca^{2+}\) sensitivity of secretion (see Table 1) and \(Ca_{null}\) represents the minimal \(Ca^{2+}\) concentration necessary for insulin release.

For one dimensional case, the equation (1) is given by:

$$\begin{aligned} \begin{aligned} \frac{\partial ^{2}[Ca^{2+}]}{\partial x^2}-\frac{k_{i} ^{+}[B_{i}]_{\infty }}{D_{Ca}}([Ca^{2+}]-[Ca^{2+}]_{\infty })+\frac{P_{ER}}{D_{Ca}}([Ca^{2+}]_{ER}-[Ca^{2+}])\\frac{P _{SERCA}^{max}}{D_{Ca}}\frac{[Ca^{2+}]^2}{k_{pump} ^{2}+[Ca^{2+}]^2}=\frac{1}{D_{Ca}}\frac{\partial [Ca^{2+}]}{\partial t}. \end{aligned} \end{aligned}$$
(10)

2.1 Initial and Bounadry Conditions

The \(\beta\)-cells are round and oval in shape with an average diameter 10 \(\mu m\) (Zimny and Blackard 1975). A linear cell is assumed to have a length of x=10 \(\mu m\) for a one-dimensional single-cell model. It is assumed that the calcium source is present at x=0 \(\mu m\) and therefore, flux boundary condition for the source is given by Smith et al. (1996); Smith (1996):

$$\begin{aligned} \lim _{x\rightarrow 0}\left( -D_{Ca}\frac{\partial [Ca^{2+}]}{\partial x}\right) =\sigma _{Ca}. \end{aligned}$$
(11)

The other end of the boundary is maintained at the background (resting) cytosolic calcium concentration (Bertram et al. 2006; Felix-Martinez and Godinez-Fernandez 2017) and expressed as given below:

$$\begin{aligned} \lim _{x\rightarrow 10} [Ca^{2+}]=0.1 \mu M. \end{aligned}$$
(12)

At the initial time t = 0 msec, the system is supposed to be at rest and thus, the initial calcium concentration at rest is given by:

$$\begin{aligned} \lim _{t\rightarrow 0}[Ca^{2+}]=0.1 \mu M. \end{aligned}$$
(13)

The study can be performed for two cases since the non-linear terms are negligible (Fig. 1). Therefore, the terms of SERCA pump can be linearised in the form of the following two cases:

Fig. 1
figure 1

Variation in the linear and non-linear behavior of SERCA pump flux along with the calcium concentration

Full size image

  • CASE-I When       \(k_{pump}>>[Ca^{2+}]\), then

    $$\begin{aligned} \frac{[Ca^{2+}]^2}{k_{pump} ^{2}+[Ca^{2+}]^2}\le \frac{[Ca^{2+}]^2}{k_{pump} ^{2}}\le \frac{[Ca^{2+}] }{k_{pump} }. \end{aligned}$$
    (14)

Equation (10) can be written as:

$$\begin{aligned} \begin{aligned} \frac{\partial ^{2}[Ca^{2+}]}{\partial x^2}-\frac{k_{i} ^{+}[B_{i}]_{\infty }}{D_{Ca}}([Ca^{2+}]-[Ca^{2+}]_{\infty })+\frac{P_{ER}}{D_{Ca}}([Ca^{2+}]_{ER}-[Ca^{2+}])\\ -\frac{P _{SERCA}^{max}}{D_{Ca}}\frac{[Ca^{2+}] }{k_{pump} }=\frac{1}{D_{Ca}}\frac{\partial [Ca^{2+}]}{\partial t}. \end{aligned} \end{aligned}$$
(15)

Rearranging the terms:

$$\begin{aligned} \begin{aligned} \frac{\partial ^{2}[Ca^{2+}]}{\partial x^2}-\frac{1}{D_{Ca}}\left( k _{i}^{+}[B_{i}]_{\infty }+P_{ER}+\frac{P _{SERCA}^{max}}{k_{pump}}\right) [Ca^{2+}]\\ +\frac{1}{D_{Ca}}\left( k_{i}^{+}[B_{i}]_{\infty }[Ca^{2+}]_{\infty }+P_{ER}[Ca^{2+}]_{ER}\right) =\frac{1}{D_{Ca}}\frac{\partial [Ca^{2+}]}{\partial t}, \end{aligned} \end{aligned}$$
$$\begin{aligned} \frac{\partial ^{2}[Ca^{2+}]}{\partial x^2}-a_1[Ca^{2+}]+b_1-\frac{1}{D_{Ca}}\frac{\partial [Ca^{2+}]}{\partial t}=0, \end{aligned}$$
(16)

where

$$\begin{aligned} a_1=\frac{1}{D_{Ca}}\left( k_{i} ^{+}[B_{i}]_{\infty }+P_{ER}+\frac{P _{SERCA}^{max}}{k_{pump}}\right) , \end{aligned}$$
(17)

and

$$\begin{aligned} b_1= \frac{1}{D_{Ca}}\left( k_{i}^{+}[B_{i}]_{\infty }[Ca^{2+}]_{\infty }+P_{ER}[Ca^{2+}]_{ER}\right) . \end{aligned}$$
(18)

  • CASE-II When       \(k_{pump}<<[Ca^{2+}]\), then       \(k_{pump}=\lambda [Ca^{2+}]\)      for      \(0<\lambda <1\),

    $$\begin{aligned} \frac{[Ca^{2+}]^2}{(k_{pump} ^{2}+ [Ca^{2+}]^2)}=\frac{[Ca^{2+}]^2}{(\lambda ^{2}[Ca^{2+}] ^{2}+[Ca^{2+}] ^{2})}= \frac{1}{\lambda ^{2}+1 }. \end{aligned}$$
    (19)

Substituting equation (19) in equation (10) we get,

$$\begin{aligned} \begin{aligned} \frac{\partial ^{2}[Ca^{2+}]}{\partial x^2}-\frac{k_{i}^{+}[B_{i}]_{\infty }}{D_{Ca}}([Ca^{2+}]-[Ca^{2+}]_{\infty })+\frac{P_{ER}}{D_{Ca}}([Ca^{2+}]_{ER}-[Ca^{2+}])\\ -\frac{P _{SERCA}^{max}}{D_{Ca}}\frac{1}{\lambda ^{2}+1} =\frac{1}{D_{Ca}}\frac{\partial [Ca^{2+}]}{\partial t}. \end{aligned} \end{aligned}$$
(20)

Rearranging the terms:

$$\begin{aligned}{} & {} \frac{\partial ^{2}[Ca^{2+}]}{\partial x^2}-\frac{1}{D_{Ca}}\left( k_{i}^{+}[B_{i}]_{\infty }+P_{ER}\right) [Ca^{2+}]\nonumber \\{} & {} +\frac{1}{D_{Ca}}\left( k_{i} ^{+}[B_{i}]_{\infty }[Ca^{2+}]_{\infty }+P_{ER}[Ca^{2+}]_{ER}-P _{SERCA}^{max}\frac{1}{\lambda ^{2}+1} \right) =\frac{1}{D_{Ca}}\frac{\partial [Ca^{2+}]}{\partial t}, \end{aligned}$$
(21)
$$\begin{aligned}{} & {} \frac{\partial ^{2}[Ca^{2+}]}{\partial x^2}-a_2[Ca^{2+}]+b_2-\frac{1}{D_{Ca}}\frac{\partial [Ca^{2+}]}{\partial t}=0, \end{aligned}$$
(22)

where

$$\begin{aligned} a_2=\frac{1}{D_{Ca}}\left( k_{i}^{+}[B_{i}]_{\infty }+P_{ER}\right) , \end{aligned}$$
(23)

and

$$\begin{aligned} b_2= \frac{1}{D_{Ca}}\left( k_{i}^{+}[B_{i}]_{\infty }[Ca^{2+}]_{\infty }+P_{ER}[Ca^{2+}]_{ER}-P _{SERCA}^{max}\frac{1}{\lambda ^{2}+1}\right) . \end{aligned}$$
(24)

From equations (16) and (22), we can write:

$$\begin{aligned} \frac{\partial ^{2}[Ca^{2+}]}{\partial x^2}-a[Ca^{2+}]+b-\frac{1}{D_{Ca}}\frac{\partial [Ca^{2+}]}{\partial t}=0, \end{aligned}$$
(25)

where

      for Case-I: \(a=a_1\) and \(b=b_1\), and

      for Case-II: \(a=a_2\) and \(b=b_2\).

The numerical solution of equation (25) was obtained using the finite element method. The cell is discretized into 40 elements, and a piecewise linear shape function is used for interpolation (See Appendix A for the solution procedure). The Crank-Nicolson scheme is used for handling the time derivative terms and simulated on MATLAB to obtain numerical results. The Numerical data of various parameters employed for simulation are given in Table 1.

Table 1 Biophysical parameters and numerical data
Full size table

3 Results and Discussion

Numerical results have been computed using the data of parameters given in Table 1. The profiles for calcium concentration, \(IP_{3}\) production, degradation and net growth of \(IP_{3}\), production of ATP and insulin secretion have been plotted.

Figure 2 displays the spatio-temporal \(Ca^{2+}\) distribution in a \(\beta\)-cell for standard values listed in Table 1. Figure 2A represents the \(Ca^{2+}\) distribution along spatial dimensions at different points of time t = 0 msec, 20 msec, 35 msec, 70 msec, 145 msec and 5000 msec, respectively. In Fig. 2A \(Ca^{2+}\) concentration is initially high at the source and falls along spatial dimensions and attains the value 0.1 \(\mu M\) on the other end. This fall in concentration profiles occurs due to SERCA pump activity, cytosolic buffers and diffusion coefficients. The calcium profiles show non-linear behavior along the space dimension, and these curves approach the linear behavior with an increase in time. This is because the dynamic processes of calcium signaling are trying to adjust with each other for its regulatory functions. With an increase in time, the processes approach equilibrium, making the curves close to linear behavior.

Fig. 2
figure 2

Calcium dynamics for standard values from Table 1

Full size image

In Fig. 2B, the calcium concentration at different nodal values x = 0, 2.50, 5.0, 7.50 and 10.0 \(\mu m\) increases with time and attains a steady state after 500 msec. This increase in calcium concentration occurs due to the activity of source influx and ER-leak. The calcium level is higher at x = 0 \(\mu m\) and lowest at x=10.0 \(\mu m\). This is due to the activity of the SERCA pump, buffer and diffusion coefficient. The profiles show non-linear behavior along the temporal dimension.

Figure 3 shows the spatio-temporal calcium-dependent \(IP_{3}\)-production and degradation in a \(\beta\)-cell. Figure 3A and Fig. 3B represent the \(IP_{3}\) production at different points of time t = 0 msec, 20 msec, 35 msec, 70 msec, 145 msec and 5000 msec, and different positions x = 0, 2.50, 5.0, 7.50 and 10.0 \(\mu m\), respectively. It can be observed from the figure that \(IP_{3}\)-production shows behavior similar to the calcium concentration profile. The \(IP_{3}\)-production increases with time for different spatial positions in almost the same way as seen in Fig. 2 for calcium concentration profiles. This is due to the dependence of \(IP_{3}\)-production on calcium concentration in a \(\beta\)-cell.

Fig. 3
figure 3

\(IP_{3}\)-production and degradation for standard values as given in Table 1

Full size image

Figure 3C and Fig. 3D represent the \(IP_{3}\)-degradation at different points of time t = 0 msec, 20 msec, 35 msec, 70 msec, 145 msec and 5000 msec and different positions x = 0, 2.50, 5.0, 7.50 and 10.0 \(\mu m\), respectively. In Fig. 3C, the shape of the \(IP_{3}\)-degradation curve is concave at times t=20 msec and 35 msec, which is similar to the shape of the calcium concentration profile in Fig. 2A. This concavity of the \(IP_{3}\)-degradation curve decreases with time and starts changing into a convex shape from time t = 70 msec. The convexity of the curve increases with time, which is dissimilar to the behavior of the calcium concentration profile in Fig. 2A. This departure in the behavior of \(IP_{3}\)-degradation from the behavior of calcium profiles with increased time may be because the non-linear term in equation (6) increases with calcium concentration. At a low calcium concentration level, the decrease in \(IP_{3}\)-degradation along a spatial dimension is sharp initially. Still, as the calcium concentration rises, the fall in \(IP_{3}\)-degradation along a spatial dimension near x = 0 \(\mu m\) becomes gradual and smooth, and the curve becomes convex. It is observed from Fig. 3B and Fig. 3D that the temporal curves show almost similar behavior as observed for calcium temporal curves in Fig. 2B. The main difference is the time step at which they attain their steady state. The production rate of \(IP_{3}\) attains the steady state at 300 msec, whereas the steady state for degradation of \(IP_{3}\) is attained at 150 msec for the initial position and 320 msec for the position x = 7.50 \(\mu m\).

Figure 4 exhibits the spatio-temporal calcium-dependent and \(IP_{3}\)-dependent ATP-production in a \(\beta\)-cell. Figure 4A and B represent the \(Ca^{2+}\)-dependent ATP-production at different points of time t = 0 msec, 20 msec, 35 msec, 70 msec, 145 msec and 5000 msec and different positions x = 0, 2.50, 5.0, 7.50 and 10.0 \(\mu m\), respectively. The figure shows that the ATP-production is high at the initial node and decreases gradually as the distance increases. The \(Ca^{2+}\)-dependent ATP-production in Fig. 4A at times t = 20 msec, 35 msec and 70 msec is similar to the concave curve seen in Fig. 2A for calcium profile. However, this concavity in \(Ca^{2+}\)-dependent ATP-production decreases with an increase in time and becomes almost linear after t = 145 msec. This is due to the value of the non-linear term of the \(Ca^{2+}\)-dependent ATP-production approaching the linear behavior with increased calcium concentration after some time. The temporal curves in Fig. 4B of calcium-dependent ATP are similar to those in Fig. 2B. The main difference is that the temporal curves of \(Ca^{2+}\)-dependent ATP-production achieve a steady state in time t= 300 msec, which is quite earlier than the steady state achieved by the calcium profile in Fig. 2B.

Fig. 4
figure 4

\(Ca^{2+}\) and \(IP_{3}\)-dependent ATP-production for standard values as given in Table 1

Full size image

Figure 4C and D represent the \(IP_{3}\)-dependent ATP-production at different points of time t=0 msec, 20 msec, 35 msec, 70 msec, 145 msec and 5000 msec, and different positions x = 0, 2.50, 5.0, 7.50 and 10.0 \(\mu m\) respectively. The \(IP_{3}\)-dependent ATP-production curves show the concave behavior at time t = 20 msec, 35 msec and 70 msec, which is similar to the calcium profiles and IP\(_{3}\)-production in Figs. 2A and 3A. However, the behavior of \(IP_{3}\)-dependent ATP-production changes with an increase in time and becomes convex after t = 145 msec. The fall in curves is sharp around x = 0 \(\mu m\) where the \(Ca^{2+}\) concentration is low. Still, with increased calcium concentration, the \(IP_{3}\)-dependent ATP-production falls gradually near x=0 \(\mu m\) and becomes sharper near the other end. In Fig. 4D, the temporal profile of \(IP_{3}\)-dependent ATP-production shows the same behavior as the temporal calcium concentration and \(IP_{3}\) profiles in Figs. 2B and Fig. 3B with the difference in their time of achieving steady states. The \(IP_{3}\)-dependent ATP-production steady state reaches the time of 300 msec.

Figure 5 shows the net growth of \(IP_{3}\) with respect to space and time for standard values, as listed in Table 1. Figure 5A displays the \(IP_{3}\) net growth along spacial dimensions at different points of time t = 0 msec, 20 msec, 35 msec, 70 msec, 145 msec and 5000 msec, respectively. It is observed from Fig. 5A that at the first node x = 0 \(\mu m\), the net growth is high for t= 5000 msec and low for the initial time. As we move far from the source, net growth starts decreasing and converges to the initial concentration of \(IP_{3}\) (0.16 \(\mu M\)) at the other end. Figure 5B shows the net growth of \(IP_{3}\) with respect to time for different positions x = 0, 2.50, 5.0, 7.50 and 10.0 \(\mu m\), respectively. It is observed that \(IP_3\) increases with respect to time and attains a steady state after t = 400 msec, which is almost similar to the calcium concentration behavior as observed in Fig. 2B.

Fig. 5
figure 5

Net growth of \(IP_{3}\) for standard parameters as listed in Table 1

Full size image

Figure 6 displays the spatio-temporal changes in the amount of insulin secreted by the \(\beta\)-cell for standard values listed in Table 1. Figure 6A exhibits the insulin secretion along space for different points of time t = 0 msec, 20 msec, 35 msec, 70 msec, 145 msec and 5000 msec, respectively. It can be observed from Fig. 6A that insulin secretion at x = 0 \(\mu m\) is high and decreases gradually along the space dimension. The shape of the curves is concave and similar to the calcium profile in Fig. 2A. This is due to the fact that insulin secretion depends upon calcium concentration. When \(Ca^{2+}\) concentration is more significant than 0.1 \(\mu M\) in the cell, the insulin will be continuously secreted and the secretion rate will depend upon the cell’s sensitivity. However, calcium concentration can not always remain high in the cell as high calcium is toxic to the cell. Therefore, the calcium regulation mechanism reduces the calcium level so that the insulin secretion stops. This implies that insulin secretion occurs in switch-on and switch-off modes in the cell. This is also evident from the fact that, at the same time, all the insulin granules are not in secretion mode. There are \(\sim\)10000 granules in a \(\beta\)-cell, out of which only one hundred are in active mode, secreting insulin (Bratanova-Tochkova et al. 2002; Wang and Sherman 2008). This implies that turn by turn, the \(\beta\)-cells become sensitive for secretion in batches and are in switch-on and switch-off mode for insulin secretion.

Fig. 6
figure 6

Insulin secretion for standard parameters as listed in Table 1

Full size image

Figure 6B represents the variation in insulin secretion rate with respect to time for different positions x = 0, 2.50, 5.0, 7.50 and 10.0 \(\mu m\), respectively. It is observed from the graphs that the secretion rate of insulin with respect to time is linear. The secretion rate grows faster for the first node x = 0 \(\mu m\) than the other nodes.

Figure 7A displays the temporal variation of insulin secretion rate for changes in values of the buffer. The insulin secretion rate decreases with the increase in values of the buffer. This occurs due to the binding of buffers with free calcium reducing the calcium concentration of cells which in turn brings down the insulin secretion rate. Figure 7B exhibits the temporal changes in insulin secretion for different values of source influx and standard values of other parameters. One can notice that insulin secretion increases with the increase in source influx. This occurs due to an increase in source influx, causing an elevation in calcium concentration which in turn elevates the insulin secretion rate.

Fig. 7
figure 7

Temporal variation of insulin secretion for different values of EGTA and source influx and spatio-temporal variation of insulin secretion for different values of \(Ca^{2+}\) sensitivity of secretion

Full size image

Figure 7C and D display the spatial and temporal insulin secretion curves for different values of calcium sensitivity. The insulin secretion is maximum at x = 0 \(\mu m\) in Fig. 7C and decreases along the space towards the other end of the cell. The insulin secretion rate is high for the high values of calcium sensitivity of secretion. In Fig. 7D, insulin secretion increases with increased calcium sensitivity of secretion. This implies that the calcium sensitivity of secretion causes the granules to be in secretion mode. Therefore, for higher values of calcium sensitivity of secretion, the higher number of granules will be in secretion mode, thereby increasing insulin secretion.

Figure 8 displays the spatial and temporal calcium concentration changes for \(\sigma\)=0, 1, 10, 15 and 20 pA, respectively. Figure 8A shows the variation of calcium at t = 0.5 sec for different values of source influx. It can be observed from Fig. 8A that as the source influx increases, calcium concentration also increases at the initial node x = 0 \(\mu m\) and decreases rapidly towards the endpoint, and attains an equilibrium value 0.1 \(\mu M\) at x = 10.0 \(\mu m\).

Fig. 8
figure 8

Calcium dynamics for \([B_{i}]_{\infty }=5 \mu M\) and \(\sigma\) = 0, 1, 10, 15 and 20 pA respectively

Full size image

Figure 8B displays the variation in calcium profiles along time for various values of source influx at x = 0 \(\mu m\). Graphs show that the \(Ca^{2+}\) concentration increases and reaches the steady state at t = 200 msec for \(\sigma\)=1 pA and at t = 400 msec for \(\sigma\)=20 pA. It means that for a higher value of source influx, more time is required to attain the steady state.

Figure 9 represents the \(IP_{3}\) production and degradation with respect to time and space for \(\sigma\)=0, 1, 10, 15 and 20 pA, respectively. Figure 9A shows the behavior of \(IP_{3}\)-production at time t=0.5 sec with respect to space for low and high rates of source influx. It is observed that \(IP_{3}\)-production has a similar behavior as seen in Fig. 8A for the \(Ca^{2+}\) concentration profile. The temporal behavior of \(IP_{3}\)-production for x=0 \(\mu m\) is shown in Fig. 9B, which also possesses a similar behavior of temporal calcium profile in Fig. 8B.

Fig. 9
figure 9

\(IP_{3}\) production and degradation for \([B_{i}]_{\infty }\) = 5 \(\mu M\) and \(\sigma\) = 0, 1, 10, 15 and 20 pA respectively

Full size image

Figure 9C displays the degradation rate of \(IP_{3}\) at time t = 0.5 sec. In Fig. 9C, the behavior of the shape of \(IP_{3}\) degradation curve is concave for \(\sigma\)= 0 pA and 1 pA, and graphs start changing their nature from concave to convex as the value of source influx increases. Figure 9D shows the degradation rate of \(IP_{3}\) at x = 0 \(\mu m\) with respect to time. It can be concluded from the graphs that if the source influx is high, the degradation rate increases rapidly and attains a steady state faster as compared to the low value of source influx. It can be observed from the graphs that \(IP_{3}\)-production and degradation rate is high for larger values of source influx.

Figure 10 exhibits the \(Ca^{2+}\) and \(IP_{3}\)-dependent ATP-production in a \(\beta\)-cell for \([B_{i}]_{\infty }\) = 5 \(\mu M\) and \(\sigma\) = 0, 1, 10, 15 and 20 pA, respectively. Figure 10A shows the variation in \(Ca^{2+}\)-dependent ATP-production at time t = 0.5 sec with respect to the spatial dimension. It is observed from the Fig. 10A that for the low value of sigma ATP-production is low, and with the increment in sigma value, production starts increasing and along the space converging to a fixed quantity. This phenomenon occurs since the source influx increases at the initial node, the calcium concentration increases and so does the ATP-production. Figure 10B represents the \(Ca^{2+}\)-dependent ATP-production along the time dimension for x = 0 \(\mu m\). The nature of temporal curves is similar to the curves, as shown in Fig. 8B; the only difference is that the \(Ca^{2+}\)-dependent ATP-production curves attain the steady state earlier at t = 200 msec as compared to the calcium curves.

Fig. 10
figure 10

\(Ca^{2+}\) and \(IP_{3}\)-dependent ATP-production for \([B_{i}]_{\infty }\) = 5 \(\mu M\) and \(\sigma\) = 0, 1, 10, 15 and 20 pA respectively

Full size image

Figure 10C shows the \(IP_{3}\)-dependent ATP-production with respect to space for t = 0.5 sec, and it is observed that \(IP_{3}\)-dependent ATP-production increases as the source influx increases. The curves show the linear behavior for \(\sigma\)= 1 pA and 10 pA, and for the higher value of source influx \(IP_{3}\)-dependent ATP-production curves start changing their behavior and become convex. \(IP_{3}\)-dependent ATP-production along time for x=0 \(\mu m\) is shown in Fig. 10D. It is clearly visible in the figure that with the increase in source influx, the \(IP_{3}\)-dependent ATP-production increases rapidly and attains a steady state earlier for high values of source influx and takes more time for low values of source influx.

Figure 11 displays the spatial and temporal changes in \(Ca^{2+}\) concentration for low source influx \(\sigma\) = 1 pA and high buffer concentrations \([B_{i}]_{\infty }\) = 20 \(\mu M\). It can be observed from Fig. 11A that along the space, the curves start overlapping as time increases, indicating that the calcium concentration becomes time-independent. However, the calcium concentration increases with respect to time for initial node x = 0 \(\mu m\). This occurs due to the reason that the high concentration of buffer binds all the free calcium. The influence of high buffer concentration on calcium profile with respect to time at different positions can be observed from Fig. 11B. The figure shows that for the high value of buffer concentration, calcium starts oscillating initially and attains the steady state after a fixed time t = 400 msec approximately.

Fig. 11
figure 11

Calcium dynamics for \(\sigma =1 pA\) and \([B_{i}]_{\infty }=20\mu M\)

Full size image

Figure 12 shows the temporal variation of \(Ca^{2+}\) concentration, \(IP_{3}\)-production, \(IP_{3}\)-degradation, net growth of \(IP_{3}\), \(Ca^{2+}\)-dependent ATP-production and \(IP_{3}\)-dependent ATP-production for \(\sigma =15 pA\) and \([B_i]_{\infty }=100 \mu M\), \([B_i]_{\infty }=150 \mu M\) and \([B_i]_{\infty }=200 \mu M\) respectively. Figures show that temporal concentration increases with time. It can be noticed that for large values of buffer concentrations, the oscillations in calcium concentration, \(IP_{3}\)-production, \(IP_{3}\)-degradation, net growth of \(IP_{3}\), \(Ca^{2+}\)-dependent ATP-production, and \(IP_{3}\)-dependent ATP-production initially converge to the steady state after some time. The oscillations in these profiles are larger for larger values of buffers.

Fig. 12
figure 12

Variation of concentrations of calcium profile, \(IP_{3}\)-production, \(IP_{3}\)-degradation, net growth of \(IP_3\), \(Ca^{2+}\)-dependent ATP production and \(IP_{3}\)-dependent ATP-production for \(\sigma =15 pA\) and \([B_i]_{\infty }=100 \mu M\), \([B_i]_{\infty }=150 \mu M\) and \([B_i]_{\infty }=200 \mu M\) respectively

Full size image

4 Conclusion

The spatio-temporal mechanisms of calcium regulating \(IP_{3}\), ATP and insulin production have been modeled effectively to explore relationships among these mechanisms in a pancreatic \(\beta\)-cell. The linear finite element method and Crank-Nicolson schemes are used respectively with respect to space and time coordinates. The absolute relative approximate error (Kaw et al. 2011) has been computed to estimate the accuracy of the model given in Table 2. The maximum absolute relative approximate error was found to be 0.20%, and the minimum accuracy was found to be 99.80%. The spectral radius has been calculated to find out the stability of the system (A14) and was found to be 0.9794, which is less than 1 and thus leads to the conclusion that the method of solution is stable (Özis et al. 2003). The method used in the current study is very effective and useful as it allows to solve the complex model with the minimum effort. The FEM scheme is superior to other finite difference schemes as it provides piecewise approximation where as finite differences provide only point wise approximation. FEM provides more flexibility as compared to other numerical schemes in the present conditions of the problem. With the help of FEM, it is possible to incorporate all the parameters such as buffers, leak and SERCA pump in the model and analyze their impacts effectively. No experimental results for the present condition of the problem are available for validation. However, the results agree with the biological principles of the mechanisms occurring in the cell. The results of the numerical solution are used to obtain the conclusions of the study as given below:

  1. 1.

    The model effectively provides the cell’s response time in achieving the maximal production rates or degradation rate, etc.

  2. 2.

    Source influx and buffers effectively maintain the calcium dynamics in a \(\beta\)-cell and other factors involved in the mechanisms.

  3. 3.

    The changes in ATP-production rates follow similar changes occurring in calcium concentration due to changes in source influx, buffers, etc. Similar results for ATP concentration are reported in the findings of Ainscow and Rutter (2002).

  4. 4.

    The model is able to effectively predict the changes in \(IP_{3}\)-production rates, degradation rates, and net growth rate due to calcium concentration. It can be concluded from the model that the behavior of \(IP_{3}\) profile depends upon the calcium profiles, which is similar to the findings of Politi et al. in the sense that the alteration in calcium concentration causes alterations in \(IP_{3}\) concentration. In other words, the change in \(IP_{3}\) concentration is proportional to or similar to the change in calcium concentration in the cell in our work as well as in the work by Politi et al. (2006).

  5. 5.

    The changes in insulin secretion rate follow similar changes in \(Ca^{2+}\) concentration caused by the changes in source influx and buffers, etc. The level of insulin is found to be low in Type-1 and Type-2 diabetes (Jaberi-Douraki et al. 2015; Rorsman 2005). The insulin secretion is reported in the range of 1600 pmol/min to 400 pmol/min in normal conditions (Rorsman 2005). From Fig. 6, it can be noticed that the values of insulin secretion lie within this range. As the source influx increases, the insulin secretion elevates, which can cause the over-stimulation of the \(\beta\)-cell and lead to obesity, cardiovascular diseases, hypertension and other metabolic disorders (Sprietsma and Schuitemaker 1994).

  6. 6.

    The model is applicable in the study of Type-1 and Type-2 diabetes, as it gives information about the factors involved in insulin secretion. It is clear from the model that decreases in source influx, increase in buffer activity, diffusion coefficient and pumping rate of SERCA can decrease calcium concentration and insulin secretion below the required level in a \(\beta\)-cell, which can cause Type-2 diabetes. Some studies are reported that claim insufficient levels of insulin can cause Type-2 diabetes. On the other hand, Type-1 diabetes is an autoimmune disease characterized by the progressive destruction of \(\beta\)-cells. From the model, it is clear that dysregulation in any factor like buffer, calcium, diffusion coefficient, ER leak and SERCA pump can cause a disturbance in the whole signaling process, which can initiate \(\beta\)-cell death and can cause Type-1 diabetes.

Thus, it can be concluded that dysfunction of calcium source influx and buffers can cause high and low levels of insulin in the cell, leading to metabolic disorders, diabetes, obesity, hypertension, etc. The novel findings of the present work are as follows:

  1. 1.

    The results provide information about the time frame required by the cell to achieve the specific levels of production of \(IP_{3}\), ATP and insulin, as well as degradation of \(IP_{3}\) in the \(\beta\)-cell.

  2. 2.

    The alterations in source influx and buffers can cause changes in calcium concentration levels, production of \(IP_{3}\), ATP and secretion of insulin and \(IP_{3}\)-degradation in the \(\beta\)-cell.

The proposed model provides crucial information about the range of changes in processes of calcium dynamics causing a proportionate range of changes in \(IP_3\), ATP and insulin in a pancreatic \(\beta\)-cell. This information can be used to estimate the proportionate changes in the functions of the cell. In all, the model gives us interesting information that will be useful for optimizing the design of protocols and devices for understanding disease progression, devising techniques for diagnosis, therapies as well as predicting the response to treatments. Finally, it can be concluded that the failure or dysfunction of the mechanisms of \(Ca^{2+}\) regulation can be the cause of an increase or decrease in insulin secretion rate, ATP-production rate, \(IP_{3}\)-production, and degradation rate leading to various metabolic disorders, diabetes, obesity, hypertension, etc.

Table 2 Error Analysis of calcium concentration at x=0 \(\mu M\) for 30 and 40 elements
Full size table