|
nif:isString
|
-
We have considered cylindrical single cell morphology for our model where the cylinder length and diameter values (Table 1) are chosen to represent a DSM cell. Table 1 also provides membrane capacitance (Cm), membrane resistance (Rm) and axial resistance values for our model. The individual membrane current components that were modeled were (i) three inward currents: L-type and T-type Ca2+ currents (ICaL and ICaT), a hyperpolarization-activated current (Ih); (ii) six outward K+ currents: two voltage-gated K+ currents (IKv1, IKCNQ), an ATP-dependent K+ current (IKATP) and three Ca2+-activated K+ currents (IBK(Ca), ISK(Ca), IIK(Ca)); and (iii) an outward back-ground leak current (ILeak).
Table data removed from full text. Table identifier and caption: 10.1371/journal.pone.0200712.t001 Values of model parameters. To generate a calcium transient, this model also incorporates a simple calcium dynamics based on exponential function (details in last section of method). Formulation of a conceptual model, which is expressed in a mathematical form, is the first step in translation of a physical system to a computational model. Here the conceptual model for individual ionic current is based on the classical Hodgkin-Huxley approach [47]. The cell membrane is described as an equivalent electrical circuit consisting of a membrane capacitance connected in parallel with a number of variable conductances representing the ion channels. Fig 1 presents a parallel conductance model consisting of a Cm shunted by a variety of ion channel conductances gion with respective Nernst potential Eion.
Figure data removed from full text. Figure identifier and caption: 10.1371/journal.pone.0200712.g001 A DSM cell parallel conductance model.It consists of voltage gated Ca2+ channels, Voltage gated K+ channels, Ca2+ activated K+ channels and leakage currents. Applying Kirchhoff’s current law after injecting stimulus current Istim, we get the following differential equation describing changes in transmembrane potential Vm.
dVmdt=−1Cm(ICa+IK+Ih+Ileak+Istim)(1)All membrane currents except large conductance (BK) Ca2+-dependent activation K+ channel were modeled using the Hodgkin-Huxley formalism: I=g¯[m(Vm,t,[Ca2+]i)]xh[(Vm,t,[Ca2+]i)]y(Vm−Erev)(2) where g¯ is maximum ionic conductance, Erev is the ion’s reversal potential, the dimensionless gating variable ‘m’ describes the time/voltage/Ca2+-dependent activation and ‘h’ is the time/voltage/Ca2+-dependent inactivation of the channel conductance. The ‘x’ and ‘y’ are power to the functions. The variation of each gating variable (m or h) can be expressed by first order differential Eqs (3 and 4) dm(Vm,t)dt=m∞(Vm)−m(Vm,t)τm(3) dh(Vm,t)dt=h∞(Vm)−h(Vm,t)τh(4) where m∞ and h∞ are the steady-state values, τm and τh the time constants, all being functions of voltage and/or intracellular Ca2+ ionic concentrations. Here the state parameter dependence on vm for ion channels is described by the Boltzman equation m∞(Vm,t)=11+exp((Vm+Vm12)/Sm)(5) h∞(Vm,t)=11+exp((Vm+Vh12)/Sh)(6) Where V1/2 is the half activation potential and S is the slope factor. BK channels kinetics have been described by a 10-state Markov model (MM) according to a model developed by Cox et al., 1997 [48] and Cox 2014 [49] in which the channel’s Ca2+-dependence is modelled at a finer grain, thus affording greater accuracy on this front than the HH formalism. We therefore adapted this multi-state Markov model for the BK conductance. In this model there are five closed “horizontal” conformation states, namely C0, C1, C2, C3 and C4. Similarly, there are five open-oriented “horizontal” conformation states O0, O1, O2, O3 and O4, each corresponding to the appropriate closed state. The MM topology description includes the cooperative Ca2+ binding among the states to illustrate the interaction of four alpha subunits of the BK channel. The horizontal transitions among the closed states (C0, C1, C2, C3, C4) and among the open states (O0, O1, O2, O3, O4) are Ca2+ dependent. In comparison, vertical transitions occur between corresponding pairs of closed and open states (i.e., C0—O0, C1—O1, C2—O2, C3—O3) and are voltage dependent. A special case is constituted by the vertical transitions between C4 and O4, which are both voltage and Ca2+ dependent. All parameters and their values are given in the supplementary document. The open conformation state O4 permits the flow of K+ ions through the BK channels under the instantaneous electrochemical driving force (EDF). The BK current, IBK is calculated by the following equation IBK=gBK¯*O*(V−EK)(7) where gBK¯ is the maximum conductance and O is summation of O1, O2, O3 and O4.
In order to describe the calcium-dependent gating of Ca2+-dependent potassium channels and to update the equilibrium potential of the Ca2+ ion, it was necessary to calculate the intracellular Ca2+ concentration. Since the parameters governing many important factors required in order to describe intracellular Ca2+ handling are not known, including its diffusion, buffering, release, Na+-Ca2+ exchanger, and pump extrusion, we did not incorporate a biophysically detailed realistic intracellular Ca2+ dynamics in our model. Instead, we assumed that the Ca2+ which enters via Ca2+ channels instantaneously diffused within a thin sub-membrane shell and that determining the decay of [Ca2+]i could be lumped into a single-exponential function [50]. d[Ca2+]idt=−(1000*iCa2*F*d)−([Ca2+]i∞−[Ca2+]iτr)(8) where, iCa is the inward Ca2+ flux due to voltage gated Ca2+ channels, d is the depth of the sub-membrane shell, [Ca2+]i∞ is the baseline Ca2+ concentration, F is the Faraday’s constant and τr is the time constant.
The simplest model used for synaptic input assumes an instantaneous rise of the synaptic conductance gsyn(t) from 0 to maximum conductance g¯syn at time instant t0 followed by an exponential decay with a time constant τ: gsyn(t)=g¯syne−(t−t0)τ(9) Another popular synaptic model, the “alpha function” [25], describes a conductance that has a rising phase with finite rise time gsyn(t)=g¯synt−t0τe(1−(t−t0)τ)(10) However, due to a single time constant, τ, the time courses of the rise and decay are correlated and cannot be set independently. So for a physiologically realistic model, we have used a more general function describing synaptic conductance profiles consisting of a sum of two exponentials, one generating the rising and one generating the decay phase [51]. It allows these time constants to be set independently such that τ rise ≠ τ decay, and for t ≥ t0 gsyn(t)=g¯synf(e−(t−t0)τdecay−e−(t−t0)τrise)(11) The normalization factor f is included to ensure that the amplitude equals g¯syn. f=1−e−(t−t0)τrise+e−(t−t0)τdecay(12)
One-dimensional strand model for spike propagation: As outlined in the Introduction, to further test the robustness of our AP model, we ascertained whether our computational AP, validated at the single-cell level, would successfully propagate in a cable-like structure. In order to accomplish this, we started by setting up the action potential in a single cell elongated to a large length (22.2 mm). We divided the elongated cell into 111 interconnected compartments to behave as a continuous cable [42, 25], where each compartment was spatially isopotential (see Results III). The synaptic stimulus was injected at the midpoint of the cell, x = 11.1 mm and electrical activity was recorded at the point of stimulation, i.e., 11.1 mm (designated R0), and at various distances from the point of stimulation (R2), in order to characterize AP propagation. In syncytial tissues such as smooth muscle, intercellular gap junctions subserve cell-to-cell electrical communication [22, 52]. We therefore extended our model to investigate the effect of gap junction properties on propagated APs in DSM cells. Towards this end we first built a 3-cell model of electrically connected cells, incorporating a gap junction resistance, rj, between adjacent cells (see Results III), the resistance rj allowing passage of localized currents by means of point processes mechanisms.
Action potentials were induced in our DSM cell model by applying either an external stimulus current (Ist) or a current based on synaptic input (Isyn). External stimulus current was applied either as a brief rectangular pulse for single AP or with a long rectangular pulses for a series of APs. As the voltage clamp method eliminates the capacitive current, ionic currents can be studied separately. All the equations, symbols and constant parameters are defined in Supporting Information S1 Appendix, S1 Table and S2 Table. Simulations were computed using a fixed time step of 0.02 ms, using Euler Method, in a PC with an Intel (R) Core (TM) i5 CPU with 3.20 GHz dual core processor. The simulation environment used for this model is NEURON [25] used widely for realistic modelling of excitable cells. The simulation environment used for this model is NEURON [25], employed widely for realistic modelling of excitable cells at both individual and network level in computationally efficient ways. This flexible and powerful simulator creates a virtual platform to simulate a diverse range of electrophysiological activities. In NEURON, cell morphology is modelled via the use of individual sections and compartments. Membrane mechanisms are incorporated via point and distributed process. Two primary scripts in NEURON are HOC and NMODL, which enable the modelling of cell morphology (single DSM cell, long cable and 1-D network), point process mechanisms (current clamp, alpha synapse, and gap junctions), and distributed mechanisms (ion channels and calcium dynamics) respectively. Our DSM model contains a large number of parameters that must be assigned values based on the available data. Here, the majority of the parameter values have been assigned based on experimental studies. However, a limited number of free parameters, most of which are scaling factors, such as the maximum conductance values for each ionic current are modified to obtain acceptable fits to (i) ionic currents recorded under voltage clamp condition and (ii) action potentials in DSM cells. This mechanism was designed to be run at a single operating temperature 37 deg C which can be specified by the hoc assignment statement. This mechanism is also intended to be used at other temperatures, or to investigate the effects of temperature changes. The temperature sensitivity parameter “tadj” is defined as tadj=2(celsius−37)10(13) where Celsius is the "operating temperature". Having developed the model, we tested its robustness to intrinsic parameter variation. We did this by varying gmax (g¯) of each of the ionic conductances stepwise over a range of +/- 20% of its control (default) value. We observed that the simulated AP was robust to changes of this order. Thus, while AP parameters varied in the expected direction for each of the imposed changes of conductance (for instance, elevating the g¯ of the CaL resulted in elevation of AP peak amplitude, and vice versa), the AP did not “break down” under the imposition of any of these variations, i.e. it did not undergo any pathological variations in amplitude or wave shape parameters. Similar observations were obtained for variations in other intrinsic parameters such as time constants of the conductances, demonstrating the robustness of our simulated AP.
Standard error of regression (S) or root mean squared error (RMSE) is a goodness-of-fit measure we used for our fits of simulated action potentials as R2 has been found to be unsuitable for such nonlinear comparisons [53]. S is calculated by the following formula: S=Σ(YExpt−YSim)2N−K(14) where Yexpt is the experimental value, YSim is the corresponding fit value from simulation, K = number of parameters used in the fit equation (also known as the degrees of freedom), N is the number of data points. A lower value of S denotes smaller average errors and represents a good fit. We chose a value of 5% of the difference between maximum and minimum values taken by the experimental data as our threshold for a good model. S below this value is considered a good fit.
|
![[This page as RDF]](http://data.gesis.org/somesci/static/rdf-icon.gif){kind=icon}