Finite element model to study calcium distribution in oocytes involving voltage gated Ca²⁺ channel, ryanodine receptor and buffers

Peer review under responsibility of Alexandria University Faculty of Medicine.

Available online 20 March 2015

Abstract

Calcium is one of the most important signalling ions in cell biology performing numerous functions with high specificity. A calcium wave triggers life at fertilization but also can cause cell death. The means by which this single ion can be both highly specific and universal is believed to lie in its spatiotemporal dynamics mediated by ion channels, pumps, receptors and calcium buffers. During oocyte maturation the calcium signalling machinery undergoes differentiation which results in distinctly different calcium release patterns on all organizational scales from puffs to waves. The calcium concentration patterns required during different stages of oocyte maturation are still not completely known. Also the mechanisms involved in calcium dynamics in oocyte cell are still not well understood. In view of above a two dimensional model has been proposed to study calcium dynamics in an oocyte cell. The parameters such as buffers, ryanodine receptor and voltage gated calcium channel are incorporated in the model. Based on the biophysical conditions the initial and boundary conditions have been framed. The model is transformed into variational form and Ritz finite element method has been employed to obtain the solution. A program has been developed in MATLAB 7.10 for the entire problem and executed to obtain numerical results. The numerical results have been used to study the effect of buffers, RyR and VGCC on calcium distribution in oocyte. The results indicate that buffers can significantly decrease the calcium concentration and RyR & VGCC can significantly raise the calcium concentration level in the oocyte cell in order to initiate, sustain and terminate specific activities in the cell. The information generated from the model can be useful to biomedical scientists for clinical and biomedical applications.

Keywords:

• FEM
• VGCC
• Buffers
• RyR
• MATLAB
• Reaction diffusion equation

1 Introduction

Ca²⁺ is a second messenger that mediates a plethora of cellular function ranging from neurotransmitter release to fertilization. Specially Ca²⁺ signalling is encoded in the spatial, temporal and amplitude features of cytoplasmic Ca²⁺ dynamics.¹ That is in the same cell Ca²⁺ signals of disparate duration, amplitude or frequency result in different cellular response. For example localized Ca²⁺ release through ryanodine receptor in vascular smooth muscle leads to vasodilation.² Whereas global sustained Ca²⁺ signals lead to vasoconstriction.³ Ca²⁺ signals achieve this specificity by differentially activating Ca²⁺ dependent efforts based on their frequency, location, duration and amplitude. At fertilization, vertebrate eggs undergo a major transition from gametogenesis with dramatic cellular alteration referred to collectively as egg activation. Ca²⁺ is the universal signal for egg activation in all sexually reproducing species studied to date from plants to humans.^4,6 The fertilization induced Ca²⁺ signal has specific spatial and temporal dynamics which is essential to activate the egg and initiate embryonic development.^4,5 This specialized Ca²⁺ signal takes the form of a single or multiple Ca²⁺ transients depending on the species.⁴ Changes in the concentration of cytosolic free calcium have been found to be responsible for the initiation and regulation of a variety of cellular functions including cellular proliferation, secretion, metabolic, adjustments and changes in gene expression.^7,8 The spatiotemporal patterns of [Ca²⁺]_c as a result of agonist stimulation are as diverse as the roles of Ca²⁺ play in different cells. The temporal pattern of [Ca²⁺]_c observed in a variety of cells includes oscillations or repetitive spiking.^7,9–11 Some cells, most notably Xenopus oocyte also exhibit interesting spatial patterns of [Ca²⁺]_c including propagating waves and target and spiral patterns.¹² Ca²⁺ waves have also been observed in myocytes, astrocytes¹³ hepatocytes¹⁰ and airways epithelial cells.¹⁴ The dynamics of Ca²⁺ is very important in cellular physiology because Ca²⁺ regulates their activity and interactions.¹⁵ Ca²⁺ waves are dependent on the diffusion of Ca²⁺ ions both within and possibly between the cells: modulating Ca²⁺ ion diffusion may predictably alter the spatial and temporal character of the Ca²⁺ wave. Zeng and co-workers²⁵ developed a mathematical model of simulation of spontaneous Ca²⁺ oscillations in astrocytes mediated by voltage gated calcium channels (VGCC). A good number of attempts have been made by scientists on study of calcium distribution in neurons cells, astrocyte cells, but very few attempts are reported in the literature on modelling of calcium distribution in oocytes. No attempt is reported in the literature for modelling calcium distribution in oocytes in the presence of VGCC. In view of above a mathematical model has been developed to study effect of buffers and ryanodine receptor over Ca²⁺ profile in oocytes in the presence of VGCC. The model has been developed for a two dimensional unsteady state case. The finite element method¹⁷ is employed to solve the problem. A computer program has been developed in MATLAB 7.10 for the whole problem and executed on Intel(R) Core™ i3 CPU, 4.00 GB RAM, 2.40 GHz processor.

2 Mathematical model and solution

Calcium kinetics in Oocytes is governed by a set of reaction–diffusion equations which can be framed assuming the following bimolecular reaction between Ca²⁺ and buffer species^18,19(1) $[{\text{Ca}}^{2+}]+[B_j]\underset{k_j^{-}}{\overset{k_j^{+}}{\leftrightarrow }}[{\text{CaB}}_j]$(1) where [${\text{Ca}}^{2+}$], [$B_j$] and [${\text{CaB}}_j$] represent the cytosolic ${\text{Ca}}^{2+}$ concentration, free buffer concentration and calcium bound buffer concentration respectively and ‘j’ is an index over buffer species, $k_j^{+}$ and $k_j^{-}$ are on and off rates for jth buffer respectively. Using Fickian diffusion, the buffer reaction diffusion system in two dimensions is expressed as^20,30,32,33(2) $\frac{\partial [{\text{Ca}}^{2+}]}{\partial t}=D_{\text{Ca}}\left(\frac{{\partial }^2[{\text{Ca}}^{2+}]}{\partial x^2}+\frac{{\partial }^2[{\text{Ca}}^{2+}]}{\partial y^2}\right)+\sum R_j+{\sigma }_{\text{VGCC}}+{\sigma }_{\text{RyR}}+\delta (x){\sigma }_{\text{Ca}}$(2) (3) $\frac{\partial [B_j]}{\partial t}=D_{B_j}\left(\frac{{\partial }^2[B_j]}{\partial x^2}+\frac{{\partial }^2[B_j]}{\partial y^2}\right)+\sum R_j$(3) (4) $\frac{\partial [{\text{CaB}}_j]}{\partial t}=D_{{\text{CaB}}_j}\left(\frac{{\partial }^2[{\text{CaB}}_j]}{\partial x^2}+\frac{{\partial }^2[{\text{CaB}}_j]}{\partial y^2}\right)-\sum R_j$(4) where reaction term $R_j$ is given by(5) $R_j=-k_j^{+}[{\text{Ca}}^{2+}][B_j]+k_j^{-}[{\text{CaB}}_j]$(5) $D_{\text{Ca}}\text{,}{D}_{B_j}\text{,}{D}_{{\text{CaB}}_j}$ are diffusion coefficients of free calcium, free buffer and ${\text{Ca}}^{2+}$ bound buffer respectively. Where ${\sigma }_{\text{VGCC}}$ is net influx of ${\text{Ca}}^{2+}$ from the voltage gated calcium channel, ${\sigma }_{\text{RyR}}$ is net influx of ${\text{Ca}}^{2+}$ from the ryanodine receptor which is assumed to be within the cell i.e., at the centre $(x=2.5\mathrm{\mu }\text{m}y=2.5\mathrm{\mu }\text{m})$ and ${\sigma }_{\text{Ca}}$ is net influx of ${\text{Ca}}^{2+}$ from the source and $\delta (x)$ is the standard Dirac delta function placed at the Ca²⁺ source. Let ${[B_j]}_T=([B_j]+[{\text{CaB}}_j])$ be the total buffer concentration of jth buffer and the diffusion coefficient of buffer is not affected by the binding of calcium i.e., $D_{B_j}=D_{{\text{CaB}}_j}\text{.}$ Then Eq. (5(5) $R_j=-k_j^{+}[{\text{Ca}}^{2+}][B_j]+k_j^{-}[{\text{CaB}}_j]$(5) ) can be written as²¹(6) $R_j=-k_j^{+}[{\text{Ca}}^{2+}][B_j]+k_j^{-}({[B_j]}_T-[B_j])$(6)

It is assumed that the buffer concentration is present in excess inside the cytosol so that the concentration of free buffer is constant in space and time, i.e. $[B_j]\widetilde{=}{[B_j]}_{\infty }$. Under this assumption Eq. (3(3) $\frac{\partial [B_j]}{\partial t}=D_{B_j}\left(\frac{{\partial }^2[B_j]}{\partial x^2}+\frac{{\partial }^2[B_j]}{\partial y^2}\right)+\sum R_j$(3) ) is approximated by¹⁹(7) $k_j^{+}[{\text{Ca}}^{2+}][B_j]=k_j^{-}({[B_j]}_T-{[B_j]}_{\infty })$(7) where ${[B_j]}_{\infty }=\frac{k_j^{-}{[B_j]}_T}{(k_j^{-}+k_j^{+}{[{\text{Ca}}^{2+}]}_{\infty })}$ is the background buffer concentration. Thus for single mobile buffer species Eq. (2(2) $\frac{\partial [{\text{Ca}}^{2+}]}{\partial t}=D_{\text{Ca}}\left(\frac{{\partial }^2[{\text{Ca}}^{2+}]}{\partial x^2}+\frac{{\partial }^2[{\text{Ca}}^{2+}]}{\partial y^2}\right)+\sum R_j+{\sigma }_{\text{VGCC}}+{\sigma }_{\text{RyR}}+\delta (x){\sigma }_{\text{Ca}}$(2) ) can be written as^18,19(8) $\frac{\partial [{\text{Ca}}^{2+}]}{\partial t}=D_{\text{Ca}}\left(\frac{{\partial }^2[{\text{Ca}}^{2+}]}{\partial x^2}+\frac{{\partial }^2[{\text{Ca}}^{2+}]}{\partial y^2}\right)-k_j^{+}{[B_j]}_{\infty }\left([{\text{Ca}}^{2+}]-{[{\text{Ca}}^{2+}]}_{\infty }\right)+{\sigma }_{\text{VGCC}}+{\sigma }_{\text{RyR}}+\delta (x){\sigma }_{\text{Ca}}$(8) where $D_{\text{Ca}}$ is the diffusion coefficient of free calcium, $\delta (x){\sigma }_{\text{Ca}}$ is the source amplitude due to the calcium channel. ${\sigma }_{\text{VGCC}}$ is the flux due to VGCC and this has been modelled using the Goldman–Hodgkin–Kartz (GHK) current equation.^20,22 We assume a single point source of ${\text{Ca}}^{2+}$, ${\sigma }_{\text{Ca}}$ at $x=0\text{,}y=0$, there are no sources for buffers and buffer concentration is in equilibrium with ${\text{Ca}}^{2+}$ far from the source and GHK equation as(9) $I_{\text{Ca}}=P_{\text{Ca}}{z}_{\text{Ca}}^2\frac{F^2{V}_m}{\mathit{RT}}\frac{{[{\text{Ca}}^{2+}]}_i-{[{\text{Ca}}^{2+}]}_o\exp \left(-z_{\text{Ca}}\frac{{\mathit{FV}}_m}{\mathit{RT}}\right)}{1-\exp \left(-z_{\text{Ca}}\frac{{\mathit{FV}}_m}{\mathit{RT}}\right)}$(9) where [Ca²⁺]_i and [Ca²⁺]_o are the intracellular and extracellular calcium concentration respectively. $P_{\text{Ca}}$ is the permeability of calcium ion, $z_{\text{Ca}}$ is the valency of calcium ion. $F$ is the Faradays constant. $V_m$ is the membrane potential. $R$ is gas constant and $T$ is absolute temperature. Eq. (9(9) $I_{\text{Ca}}=P_{\text{Ca}}{z}_{\text{Ca}}^2\frac{F^2{V}_m}{\mathit{RT}}\frac{{[{\text{Ca}}^{2+}]}_i-{[{\text{Ca}}^{2+}]}_o\exp \left(-z_{\text{Ca}}\frac{{\mathit{FV}}_m}{\mathit{RT}}\right)}{1-\exp \left(-z_{\text{Ca}}\frac{{\mathit{FV}}_m}{\mathit{RT}}\right)}$(9) ) is converted into molar/s by using the following equation(10) ${\sigma }_{\text{VGCC}}=\frac{-I_{\text{Ca}}}{z_{\text{Ca}}{\mathit{FV}}_{\text{Oocyte}}}$(10)

The negative sign in Eq. (10(10) ${\sigma }_{\text{VGCC}}=\frac{-I_{\text{Ca}}}{z_{\text{Ca}}{\mathit{FV}}_{\text{Oocyte}}}$(10) ) is taken due to the inward current by convection. GHK current equation gives the current density as a function of voltage. The GHK equation is derived from the constant field which assumes that the electric field in the membrane is constant and thus ions move in the membrane in the same way as in free solution. ${\sigma }_{\text{RyR}}$ is flux due to ryanodine receptor given as^26,31(11) ${\sigma }_{\text{RyR}}=V_{\text{RyR}}{P}_o\left({[{\text{Ca}}^{2+}]}_{\mathit{ER}}-[{\text{Ca}}^{2+}]\right)$(11) Combining Eqs. (8(8) $\frac{\partial [{\text{Ca}}^{2+}]}{\partial t}=D_{\text{Ca}}\left(\frac{{\partial }^2[{\text{Ca}}^{2+}]}{\partial x^2}+\frac{{\partial }^2[{\text{Ca}}^{2+}]}{\partial y^2}\right)-k_j^{+}{[B_j]}_{\infty }\left([{\text{Ca}}^{2+}]-{[{\text{Ca}}^{2+}]}_{\infty }\right)+{\sigma }_{\text{VGCC}}+{\sigma }_{\text{RyR}}+\delta (x){\sigma }_{\text{Ca}}$(8) )–(11(11) ${\sigma }_{\text{RyR}}=V_{\text{RyR}}{P}_o\left({[{\text{Ca}}^{2+}]}_{\mathit{ER}}-[{\text{Ca}}^{2+}]\right)$(11) ) we get proposed mathematical model as given below(12) $\frac{\partial [{\text{Ca}}^{2+}]}{\partial t}=D_{\text{Ca}}\left(\frac{{\partial }^2[{\text{Ca}}^{2+}]}{\partial x^2}+\frac{{\partial }^2[{\text{Ca}}^{2+}]}{\partial y^2}\right)-k_j^{+}{[B_j]}_{\infty }\left([{\text{Ca}}^{2+}]-{[{\text{Ca}}^{2+}]}_{\infty }\right)+P_{\text{Ca}}{z}_{\text{Ca}}^2\frac{F^2{V}_m}{\mathit{RT}}\frac{{[{\text{Ca}}^{2+}]}_i-{[{\text{Ca}}^{2+}]}_o\exp \left(-z_{\text{Ca}}\frac{{\mathit{FV}}_m}{\mathit{RT}}\right)}{1-\exp (-z_{\text{Ca}}\frac{{\mathit{FV}}_m}{\mathit{RT}})}+V_{\text{RyR}}{P}_o\left({[{\text{Ca}}^{2+}]}_{\mathit{ER}}-[{\text{Ca}}^{2+}]\right)+\delta (x){\sigma }_{\text{Ca}}$(12)

The point source of calcium is assumed at $x=0y=0$ and as we move away from the source, the calcium concentration achieves its background value i.e., $0.1\mathrm{\mu }\text{M}$. Thus the initial and boundary conditions for the above problem are^30,32,33

Initial Condition:(13) ${[{\text{Ca}}^{2+}]}_{t=0}=0.1\mathrm{\mu }\text{M}\forall x\text{,}y$(13) Boundary Conditions:(14) ${\displaystyle \underset{x\to 0\text{,}y\to 0}{\lim }}-D_{\text{Ca}}\frac{\partial [{\text{Ca}}^{2+}]}{\partial x}={\sigma }_{\text{Ca}}$(14) (15) ${\displaystyle \underset{x\to 5\text{,}y\to 5}{\lim }}[{\text{Ca}}^{2+}]=0.1\mathrm{\mu }\text{M}$(15)

Here $[{\text{Ca}}^{2+}]$ is the background calcium concentration, $P_{\text{Ca}}[{\text{Ca}}^{2+}]$ represents the rate of calcium efflux from the cytosol into extracellular space. ${\sigma }_{\text{Ca}}$ represents the flux due to $[{\text{Ca}}^{2+}]$ and incorporated on the boundary tends to the background concentration of 0.1 $0.1\mathrm{\mu }\text{M}$ as $x\to \infty \text{,}y\to \infty$ but the domain taken by us is not infinite one. Here we are taking the distance required for $[{\text{Ca}}^{2+}]$ to attain background concentration as $5\mathrm{\mu }\text{m}$ for Oocyte along x-axis and y-axis.^30,32,33 Our problem is to solve Eq. (12(12) $\frac{\partial [{\text{Ca}}^{2+}]}{\partial t}=D_{\text{Ca}}\left(\frac{{\partial }^2[{\text{Ca}}^{2+}]}{\partial x^2}+\frac{{\partial }^2[{\text{Ca}}^{2+}]}{\partial y^2}\right)-k_j^{+}{[B_j]}_{\infty }\left([{\text{Ca}}^{2+}]-{[{\text{Ca}}^{2+}]}_{\infty }\right)+P_{\text{Ca}}{z}_{\text{Ca}}^2\frac{F^2{V}_m}{\mathit{RT}}\frac{{[{\text{Ca}}^{2+}]}_i-{[{\text{Ca}}^{2+}]}_o\exp \left(-z_{\text{Ca}}\frac{{\mathit{FV}}_m}{\mathit{RT}}\right)}{1-\exp (-z_{\text{Ca}}\frac{{\mathit{FV}}_m}{\mathit{RT}})}+V_{\text{RyR}}{P}_o\left({[{\text{Ca}}^{2+}]}_{\mathit{ER}}-[{\text{Ca}}^{2+}]\right)+\delta (x){\sigma }_{\text{Ca}}$(12) ) coupled with Eqs. (13(13) ${[{\text{Ca}}^{2+}]}_{t=0}=0.1\mathrm{\mu }\text{M}\forall x\text{,}y$(13) )–(15(15) ${\displaystyle \underset{x\to 5\text{,}y\to 5}{\lim }}[{\text{Ca}}^{2+}]=0.1\mathrm{\mu }\text{M}$(15) ). For our convenience we are writing ‘u’ in lieu of [Ca²⁺]. From (12(12) $\frac{\partial [{\text{Ca}}^{2+}]}{\partial t}=D_{\text{Ca}}\left(\frac{{\partial }^2[{\text{Ca}}^{2+}]}{\partial x^2}+\frac{{\partial }^2[{\text{Ca}}^{2+}]}{\partial y^2}\right)-k_j^{+}{[B_j]}_{\infty }\left([{\text{Ca}}^{2+}]-{[{\text{Ca}}^{2+}]}_{\infty }\right)+P_{\text{Ca}}{z}_{\text{Ca}}^2\frac{F^2{V}_m}{\mathit{RT}}\frac{{[{\text{Ca}}^{2+}]}_i-{[{\text{Ca}}^{2+}]}_o\exp \left(-z_{\text{Ca}}\frac{{\mathit{FV}}_m}{\mathit{RT}}\right)}{1-\exp (-z_{\text{Ca}}\frac{{\mathit{FV}}_m}{\mathit{RT}})}+V_{\text{RyR}}{P}_o\left({[{\text{Ca}}^{2+}]}_{\mathit{ER}}-[{\text{Ca}}^{2+}]\right)+\delta (x){\sigma }_{\text{Ca}}$(12) ) we get(16) $\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}-\alpha u+\beta -\frac{1}{D_{\text{Ca}}}\frac{\partial u}{\partial t}=0$(16) (17) $\alpha =\frac{1}{D_{\text{Ca}}}\left[k_j^{+}{[B_j]}_{\infty }-P_{\text{Ca}}{z}_{\text{Ca}}\frac{\frac{F^2{V}_m}{\mathit{RT}}}{1-\exp \left(-z_{\text{Ca}}\frac{{\mathit{FV}}_m}{\mathit{RT}}\right)}+V_{\text{RyR}}{P}_o\right]\text{,}\beta =\frac{1}{D_{\text{Ca}}}\left[k_j^{+}{[B_j]}_{\infty }{u}_{\infty }-\frac{P_{\text{Ca}}{z}_{\text{Ca}}\frac{F^2{V}_m}{\mathit{RT}}\exp \left(-z_{\text{Ca}}\frac{{\mathit{FV}}_m}{\mathit{RT}}\right)}{1-\exp \left(-z_{\text{Ca}}\frac{{\mathit{FV}}_m}{\mathit{RT}}\right)}{u}_o+V_{\text{RyR}}{P}_o{u}_{\mathit{ER}}\right]$(17)

Applying finite element method on Eq. (16(16) $\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}-\alpha u+\beta -\frac{1}{D_{\text{Ca}}}\frac{\partial u}{\partial t}=0$(16) ) we can get variational form as(18) $I^{(e)}=\frac{1}{2}\int {\int }_A\left[{\left(\frac{\partial u^{(e)}}{\partial x}\right)}^2{\left(\frac{\partial u^{(e)}}{\partial y}\right)}^2+\alpha u^{(e)2}-2\beta u^{(e)}+\frac{1}{D_{\text{Ca}}}\left(\frac{\partial {(u^{(e)})}^2}{\partial t}\right)\right]\mathit{dA}-{\mu }^{(e)}{\int }_{\partial A}\left(\frac{\sigma }{2D_{\text{Ca}}}{u}^{(e)}{|}_{x=5}\right)\mathit{dy}$(18)

Here we have used ‘u’ in lieu of [Ca²⁺] for our convenience and e = 1, 2, 3…162. Also ${\mu }^{(e)}=1$ for $e=1$ and ${\mu }^{(e)}=0$ for rest of elements. The following linear shape function for calcium concentration within each element has been taken as¹⁷(19) $u^{(e)}=c_1^{(e)}+c_2^{(e)}x+c_3^{(e)}y$(19) (20) $u^{(e)}=p^T{c}^{(e)}$(20) where(21) $\begin{array}{ccc}\hfill p^T=[1\hfill & \hfill x\hfill & \hfill y]\hfill \end{array}$(21) and(22) $\begin{array}{ccc}\hfill c^{(e)T}=[c_1^{(e)}\hfill & \hfill c_2^{(e)}\hfill & \hfill c_3^{(e)}]\hfill \end{array}$(22)

Substituting nodal conditions in Eq. (21(21) $\begin{array}{ccc}\hfill p^T=[1\hfill & \hfill x\hfill & \hfill y]\hfill \end{array}$(21) ), we get(23) ${\overline{u}}^{(e)}=P^{(e)}\ast c^{(e)}$(23) where(24) ${\overline{u}}^{(e)}=\left[\begin{array}{c}\hfill u_i\hfill \\ \hfill u_j\hfill \\ \hfill u_k\hfill \\ \hfill \hfill \end{array}\right]\text{and}{P}^{(e)}=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill x_i\hfill & \hfill y_i\hfill \\ \hfill 1\hfill & \hfill x_j\hfill & \hfill y_j\hfill \\ \hfill 1\hfill & \hfill x_k\hfill & \hfill y_k\hfill \end{array}\right]$(24)

From the Eq. (23(23) ${\overline{u}}^{(e)}=P^{(e)}\ast c^{(e)}$(23) ), we have(25) $c^{(e)}=R^{(e)}\ast {\overline{u}}^{(e)}$(25) where(26) $R^{(e)}=P^{{(e)}^{-1}}$(26)

Substituting $c^{(e)}$ from Eq. (25(25) $c^{(e)}=R^{(e)}\ast {\overline{u}}^{(e)}$(25) ) in (20(20) $u^{(e)}=p^T{c}^{(e)}$(20) ), we get(27) $u^{(e)}=p^T{R}^{(e)}{\overline{u}}^{(e)}$(27)

The expression (18(18) $I^{(e)}=\frac{1}{2}\int {\int }_A\left[{\left(\frac{\partial u^{(e)}}{\partial x}\right)}^2{\left(\frac{\partial u^{(e)}}{\partial y}\right)}^2+\alpha u^{(e)2}-2\beta u^{(e)}+\frac{1}{D_{\text{Ca}}}\left(\frac{\partial {(u^{(e)})}^2}{\partial t}\right)\right]\mathit{dA}-{\mu }^{(e)}{\int }_{\partial A}\left(\frac{\sigma }{2D_{\text{Ca}}}{u}^{(e)}{|}_{x=5}\right)\mathit{dy}$(18) ) can be written as(28) $I^{(e)}=I_l^{(e)}+I_m^{(e)}-I_n^{(e)}-I_k^{(e)}+I_p^{(e)}$(28) where(29) $I_l^{(e)}=\frac{1}{2}\int {\int }_A\left[{\left(\frac{\partial u^{(e)}}{\partial x}\right)}^2+{\left(\frac{\partial u^{(e)}}{\partial y}\right)}^2\right]\mathit{dA}$(29) (30) $I_m^{(e)}=\frac{1}{2}\alpha \int {\int }_A{u}^{{(e)}^2}\mathit{dA}$(30) (31) $I_n^{(e)}=2\beta \int {\int }_A{u}^{(e)}{u}_{\infty }\mathit{dA}$(31) (32) $I_k^{(e)}=\frac{{\mu }^{(e)}}{2}{\int }_{\partial A}\left[\frac{{\sigma }_{\text{Ca}}}{D_{\text{Ca}}}{u}^{(e)}{|}_{x=5}\right]\mathit{dy}$(32) (33) $I_p^{(e)}=\frac{1}{2D_{\text{Ca}}}\frac{d}{\mathit{dt}}\int {\int }_A{(u^{(e)})}^2\mathit{dA}$(33)

Now we extremize the integral $I^{(e)}$ w.r.t. each nodal calcium concentration $u_i$ as given below(34) $\frac{\mathit{dI}}{d\overline{u}}={\displaystyle \sum _{e=1}^{162}}{\bar{M}}^{(e)}\frac{{\mathit{dI}}^{(e)}}{d{\overline{u}}^{(e)}}{\bar{M}}^{{(e)}^T}=0$(34) and(35) $\frac{{\mathit{dI}}^{(e)}}{d{\overline{u}}^{(e)}}=\frac{{\mathit{dI}}_l^{(e)}}{d{\overline{u}}^{(e)}}+\frac{{\mathit{dI}}_m^{(e)}}{d{\overline{u}}^{(e)}}-\frac{{\mathit{dI}}_n^{(e)}}{d{\overline{u}}^{(e)}}-\frac{{\mathit{dI}}_k^{(e)}}{d{\overline{u}}^{(e)}}+\frac{d}{\mathit{dt}}\frac{{\mathit{dI}}_p^{(e)}}{d{\overline{u}}^{(e)}}$(35) where(36) ${\bar{M}}^{(e)}=\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill ·\hfill & \hfill ·\hfill & \hfill ·\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill ·\hfill & \hfill ·\hfill & \hfill ·\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill \hfill \end{array}\right]\text{and}\overline{u}=\left[\begin{array}{c}\hfill u_1\hfill \\ \hfill u_2\hfill \\ \hfill u_3\hfill \\ \hfill ·\hfill \\ \hfill ·\hfill \\ \hfill ·\hfill \\ \hfill u_{100}\hfill \end{array}\right]$(36) (37) $\text{and}I={\displaystyle \sum _{e=1}^{162}}{I}^{(e)}$(37)

This leads to a following system of linear differential equations.(38) ${[A]}_{100\times 100}{\left[\frac{\partial \overline{u}}{\partial t}\right]}_{100\times 1}+{[B]}_{100\times 100}{[\overline{u}]}_{100\times 1}={[C]}_{100\times 1}$(38) Here $[\overline{u}]=[u_1\text{,}{u}_2\text{,}{u}_3\dots u_{100}]$ A and B are the system matrices and C is the system vector. The Crank–Nicolson method is used to solve the system of differential Eq. (38(38) ${[A]}_{100\times 100}{\left[\frac{\partial \overline{u}}{\partial t}\right]}_{100\times 1}+{[B]}_{100\times 100}{[\overline{u}]}_{100\times 1}={[C]}_{100\times 1}$(38) ). A computer program has been developed in MATLAB 7.10 for the entire problem and executed on Intel(R) Core™ i3 CPU, 4.00 GB RAM, 2.40 GHz processor. The numerical values of biophysical parameters used in the model are stated in Table 1.^16,23,24,31

Table 1 Values of biophysical parameters.^16,23,24,31

Symbol	Parameter	Value
$D_{\text{Ca}}$	Diffusion coefficient	$250\mathrm{\mu }{\text{m}}^2/\text{s}$
$k_j^{+}$	On rate for EGTA	$3/\mathrm{\mu }\text{M}\text{s}$
$k_j^{-}$	Off rate for EGTA	$1/\text{s}$
$k_j^{+}$	On rate for BAPTA	$100/\mathrm{\mu }\text{M}\text{s}$
$k_j^{-}$	Off rate for BAPTA	$10/\text{s}$
${[B]}_{\infty }$	Total buffer concentration	$50\mathrm{\mu }\text{M}$
$\sigma$	Source amplitude	$1\text{pA}$
${[{\text{Ca}}^{2+}]}_o$	Extracellular Calcium Concentration	$3\text{mM}$
$P_{\text{Ca}}$	Calcium permeability	$4.3\times {10}^{-8}\text{m}/\text{s}$
$V_m$	Membrane potential	$-0.05\text{V}$
$z_{\text{Ca}}$	Valency of calcium	$2$
$V_{\text{Oocyte}}$	Volume of oocyte cytosol	$5.48\times {10}^{-11}\text{l}$
$F$	Faraday’s constant	$96487\text{C}/\text{mole}$
$R$	Gas constant	$8.314\text{J}/\text{K}\text{mole}$
$T$	Absolute temperature	$37°\text{C}$
$P_o$	Rate of calcium efflux	$0.5\text{M}/\text{s}$
$V_{\text{RyR}}$	RyR receptor rate	$0.5\mathrm{\mu }\text{M}/\text{s}$
${[{\text{Ca}}^{2+}]}_{\mathit{ER}}$	ER Ca^2+ concentration	$500\mathrm{\mu }\text{M}$

3 Results and discussion

In this part, we have shown the numerical results for calcium profile against different biophysical parameters which have been obtained using values of parameters given in Table 1 unless stated along with figures.

From Fig. 1 the effect of source amplitude and the buffers is clearly visible. We see that the calcium concentration is higher near the source and decreases slowly as we move away from the source. The concentration is higher from $0\mathrm{\mu }\text{m}\text{to}1\mathrm{\mu }\text{m}$ and then decreases slowly up to $1.5\mathrm{\mu }\text{m}$ and finally tends to its initial value of $0.1\mathrm{\mu }\text{M}$ for the source amplitude. The figure further shows that the calcium concentration is higher for lower value of buffer this is because the higher concentration of buffer binds much calcium and makes the system to attain the steady state earlier. The calcium concentration is higher from $0\mathrm{\mu }\text{m}\text{to}0.5\mathrm{\mu }\text{m}$ and then onwards decreases gradually to attain the initial value.

Figure 1 Spatial variation of calcium concentration for source amplitude $\sigma =1\text{pA}$ and for different concentration of buffer.

Fig. 2 shows the effect of parameters i.e., the calcium concentration with and without these parameters. The figure shows that the calcium concentration is lower in the absence of parameters and very much higher in the presence of them. This is clearly visible from the figure because without these parameters there is less calcium to diffuse as these parameters contribute nothing because of their absence, but the calcium concentration is higher in the presence of these parameters as they release the calcium and contribute in the diffusion so as to make the concentration high. The figure shows the Ca²⁺ concentration is low in the absence of VGCC & RyR which is $15\mathrm{\mu }\text{M}$ and increases when RyR is involved and reaches up to $25\mathrm{\mu }\text{M}$. Further the Ca²⁺ concentration reaches the higher level of $40\mathrm{\mu }\text{M}$ when VGCC is present.

Figure 2 Calcium concentration in the absence and the presence of RyR and VGCC.

Fig. 3 gives the spatial variation of calcium concentration for different time periods. We observe that at $t=50\text{ms,}100\text{ms,}150\text{ms}$ the peak calcium concentration is $27\mathrm{\mu }\text{M}28\mathrm{\mu }\text{M,}30\mathrm{\mu }\text{M}$ respectively. We observe that the peak Ca²⁺ concentration at source increases with time. Also the spread of Ca²⁺ increases with time.

Figure 3 Spatial variation of calcium concentration for different time periods.

Fig. 4 shows the steady state effect of increasing buffer concentration of EGTA buffer on calcium distribution in oocytes for source amplitude $1\text{pA}$. In all the four cases different equilibrium concentrations $({[B_j]}_{\infty })$ of EGTA buffer is taken to be $B=50\text{,}100\text{,}150\text{,}200\mathrm{\mu }\text{M}$ respectively. From the figure it is evident that as buffer concentration increases, the peak cytosolic calcium concentration at source decreases from $30\mathrm{\mu }\text{M}$ for $B=50\mathrm{\mu }\text{M}$ to $10\mathrm{\mu }\text{M}$ for $B=200\mathrm{\mu }\text{M}$. Our results suggest that if buffer concentration increases too much then it either balkanizes Ca²⁺ signal or makes it less effective. In our calculations we also observed that buffer concentration alters the time required to achieve the steady state. Here we further observe that as buffer concentration increases, the area of dispersion of cytosolic calcium concentration decreases. The results obtained in this study are in a close agreement with the experimental studies^27,28 and the results obtained by Panday and Pardasani.²⁹ and Jha et al.³⁰.

Figure 4 Spatial variation of calcium concentration for different concentrations of EGTA buffer.

4 Conclusion

The two-dimensional finite element model has been proposed for an unsteady state case and employed to study effect of VGCC, RyR and buffers on calcium concentration distribution in oocytes. The results provide us information about spatio–temporal relationship of calcium concentration with buffers, RyR and VGCC in oocyte cell. It is concluded from results that VGCC and RyR can significantly raise the level of calcium concentration in the cell in response to the requirements for initiation, sustenance of the specific activities of the cell. Also the buffers are capable of lowering down the calcium concentration in the cell when the calcium concentration becomes high in the cell for a specific activity. The cell exhibits a beautiful coordination of RyR, VGCC and buffers in order to regulate the calcium concentration in the oocyte cell. The finite element method is quite flexible and versatile in the present situation as it has been possible to incorporate the important parameters such as VGCC, RyR and buffers in the model. Such models can be developed further to study the relationship among various parameters under normal and abnormal conditions to generate the information which can be of great use to biomedical scientists for developing protocols for diagnosis and treatment of diseases related to reproduction.

Conflict of interest

The authors declare that there is no any kind of conflict of interest.

Acknowledgement

The authors are highly thankful to University Grants Commission (UGC), New Delhi, India for providing financial support to carry out this work.

Notes

Peer review under responsibility of Alexandria University Faculty of Medicine.

Available online 20 March 2015