A549 in-silico 1.0: A first computational model to simulate cell cycle dependent ion current modulation in the human lung adenocarcinoma.

Lung cancer is still a leading cause of death worldwide. In recent years, knowledge has been obtained of the mechanisms modulating ion channel kinetics and thus of cell bioelectric properties, which is promising for oncological biomarkers and targets. The complex interplay of channel expression and its consequences on malignant processes, however, is still insufficiently understood. We here introduce the first approach of an in-silico whole-cell ion current model of a cancer cell, in particular of the A549 human lung adenocarcinoma, including the main functionally expressed ion channels in the plasma membrane as so far known. This hidden Markov-based model represents the electrophysiology behind proliferation of the A549 cell, describing its rhythmic oscillation of the membrane potential able to trigger the transition between cell cycle phases, and it predicts membrane potential changes over the cell cycle provoked by targeted ion channel modulation. This first A549 in-silico cell model opens up a deeper insight and understanding of possible ion channel interactions in tumor development and progression, and is a valuable tool for simulating altered ion channel function in lung cancer electrophysiology.

Introduction

Lung cancer is one of the most prevalent forms of tumor and the leading cause of cancer death worldwide.[1-3] Increasing knowledge of molecular cancer biology and the identification of key and potentially targetable genetic and molecular aberrations that drive tumor growth provide efficient diagnostic and therapeutic approaches for lung cancer. Somatic mutations of oncogenes and tumor suppressor genes in the lung adenocarcinoma appear to be promising therapeutic oncogenic targets.[4,5] Nevertheless, despite substantial advances in early diagnosis and innovative treatment strategies survival rates still remain poor.[6] Thus, a profound understanding of cancer biology at multiple functional levels will afford novel therapeutic agents to effectively fight and cure this disease.[7] In particular, advances in our understanding of molecular alterations at genetic, epigenetic or protein expression levels together with their functional significance, and in recent years, expanding knowledge in the mechanisms and modulation of ion channel function in cancer biology lead to the development of promising cancer biomarkers and oncological targets.[8,9]

Cells are characterized by a unique composition of ion channels responsible for the bioelectric properties of the cell, playing a fundamental role in almost all cellular functions. Cancer cells, compared to their differentiated benign counterparts, typically exhibit an altered ion channel expression or activity [8,10-12] associated with tumor development and progression [10-15]. Nevertheless, no unifying pattern could as yet be identified and expression levels appear to be diverse across different cancer types.[10] The expression profile of ion channels is decisive for the membrane potential Vm, a key bioelectrical signal in the regulation of basic cellular activities such as proliferation, apoptosis, migration and differentiation.[10,14,16-18] In turn, voltage-gated ion channels (VGIC) respond to changes in the membrane potential through altered ion channel activation and conductivity, modulating the membrane potential for their part accordingly.[15,19] As a general principle, cancer cells tend to be more depolarized [19,20], whereby depolarization is assumed to initiate DNA-synthesis and mitosis.[10,16,21] During cell cycle progression, the membrane potential undergoes rhythmic oscillations starting with a further depolarization during the transition from the resting G0 to G1 phase, followed by a hyperpolarization during S phase initiation and subsequent depolarization while entering the M-G2 phase.[17,19] The exact thresholds required for driving the cells through the distinct phases have not been extensively studied and are likely to vary between different cell types.[17] It is well known, however, that the rhythmic oscillation occurs by means of a complex interplay between different predominantly hyperpolarizing, mainly voltage-dependent K+ channels (KCa, EAG, Kv, KATP K2P) and depolarizing ion channels such as the voltage-gated chloride channels (ClC) with cell cycle dependent expression levels.[16]

The crucial role of potassium conductance in governing the membrane potential and controlling the cell cycle has been supported in a number of studies.[22] It has been confirmed that blocking of potassium channels with selective inhibitors reduces proliferation in different cells.[23-29] For instance, inhibition of Kv channels reduces proliferation in prostate cancer cells [25] or analogously inhibition of Kv10.1 expression leads to reduced proliferation in diverse cancer cell lines [26]. However, it has still not been fully clarified to what extent the expression and activity of individual channels or channel-mediated changes of the membrane potential contribute to cell cycle progression, since inhibition of cell proliferation by channel blockage does not necessarily also lead to changes in the membrane potential.[10,23] Besides enhanced proliferation, expression levels of specific ion channels also correlate with other hallmarks of tumor progression such as evasion of apoptosis, sustained angiogenesis and invasion.[10,30,31] For example, Kv10.1 (EAG1) and Kv11.1 (hERG1) channels are linked to trigger angiogenesis [32-34], while inhibition of BK channels reduces migration of glioma cells [29] and the voltage-gated sodium channels Nav1.5 and Nav1.7 are generally associated with increased migration and metastasis.[10,35] By implication, targeted inhibition and activation of specific ion channels provides potential novel strategies for cancer therapy.[19,34] However, the complex interplay of channel expression, ion current dynamics together with their consequences on malignant processes are still insufficiently understood, in particular in the human lung adenocarcinoma.

In-silico modeling of whole-cell electrophysiology is a well-established tool for the description of the membrane potential in excitable cells. A range of models with different levels of complexity were introduced for simulating ion current kinetics and action potential alterations in neural or cardiac cells (e.g. neuronal model initially described by Hodgkin Huxley [36], or advanced cardiac cell models by Ten Tusscher [37] and Luo Rudy [38]). However, only a few approaches exist for non-excitable cells, e.g. the modulation of the membrane potential and cell secretion by single ion channels, calcium dynamics or activation of T-lymphocytes.[39-42]

The A549 cell line, derived from non-small cell lung cancer (NSCLC), is a widely used model for studying lung cancer and cancer drug development.[43-45] In this work we introduce for the first time an ion current model of the A549 human lung adenocarcinoma cell, representing an initial description of a cell model as a whole in cancer electrophysiology. The model takes into account the kinetics of the most relevant ion channels that contribute to the cell’s total membrane current and resting membrane potential. Based on our experimental data using the whole-cell patch-clamp technique and a comprehensive literature review (S1 Text), single channel kinetics were modeled by a hidden Markov model (HMM) approach and the number of represented ion channels estimated by fitting the macroscopic currents to the recorded whole-cell currents. The model was parameterized, taking into account the specific ion channel activities in the A549 cell obtained from the literature data and involves the main functionally expressed ion channels in the plasma membrane of the A549 cell known so far, also considering the respective voltage and calcium dependencies. This approach now allows for the first time the simulation of channel interaction, activation and inhibition, and, importantly, the prediction of membrane potential changes for parts of the cell cycle. The availability of this initial A549 in-silico model 1.0 provides a deeper understanding of the possible roles and interactions of ion channels in tumor development and progression and may aid in the testing, verification and validation of research hypotheses in lung cancer electrophysiology.

Results

Ion channels present in the A549 cell line

A review of the gene database Single Cell Expression Atlas–AMBL_EBI revealed expression of 99 ion channel genes in A549 cells. Thereunder are 52 voltage-gated ion channels, 25 ligand-gated and 22, denoted as other ion channels according to the classification of IUPHAR/BPS (guide to pharmacology) including e.g. aquaporins, sodium leak channels or the CLC family.[46] Expression of 14 additional channel genes, not listed in the AMBL_EBI, were found in an additional PubMed literature search, leading in total to 113 ion channel genes expressed in A549 cells. Fig 1 provides a summary of reported ion channels in A549 cells reviewed in this work, building the fundamental basis for the A549 in-silico model (see S1 Text for review). However, the functional expression and localization of merely 18 channels in the cell membrane of A549 was confirmed by Western Blot and patch-clamp experiments in previous studies. 11 channels were finally included in this initial cell model, comprising the fast activating and inactivating channels Kv1.3 and Kv3.4 and channels with slower, but constant activation such as TASK-1 (K2P3.1) and KCa3.1 (hIK). These channels together are decisive for the instantaneous activating current whereas Kv7.1, Kv3.1 and also KCa1.1 account for the time-dependent increase of the membrane current. The calcium entry caused by CRACM1, TRPC6 and TRPV3, and the CLC-2 mediated chloride current provoke a cumulative negative inward current, which balances the positive potassium outward current. The remaining seven to 18 channels, i.e. Kv3.3, Kv9.3, K2P9.1, ASICs, ENaC, Cav3.1 and CTFR, were not considered in the model, because the kinetics of these channels are not fully known and experimental data on them is not yet available. This level of abstraction and simplification can scarcely be regarded as representing a limitation to this “whole-cell” model, since the implemented ion channels appear to be sufficient for a valid characterization of the physiological cell function over the cell cycle.

Fig 1

Reported ion channels in the A549 cell line.

For detailed review see S1 Text.

The HMM based A549 whole-cell current model

Fig 2 illustrates the A549 whole-cell model, indicating the different ion channel types, macroscopic currents and their kinetic schemes using a hidden Markov modeling (HMM) approach. HMMs represent the gating of an ion channel through a series of conformational changes of the channel protein, assuming that the transition probability between these states depends on the present state only. Exemplarily, for the ion channel Kv7.1 a five state HMM consisting of two closed (C), two open (O) and one inactivated state (I) was implemented in accordance with the kinetic model described by Pusch et al. [47]:

Fig 2

A549 whole-cell ion current model.

The model illustrates the different ion channel types, macroscopic currents and kinetic schemes of the used hidden Markov models (HMM). Potassium channels, including Kv1.3, Kv3.1, Kv3.4, Kv7.1, TASK-1, KCA1.1 and KCa3.1 are represented in blue, green denotes the included calcium channels CRACM1, TRPV3 and TRPC6 and red represents the considered chloride channel CLC-2. See also S1 Fig. Created with BioRender.com.

The forward transition rates α, a and c, and backward transition rates β, b and d for the transitions between open and closed states (C1⇌C2⇌O1⇌O2) are voltage dependent and given in the form:

and , where α and βi represent specific gating parameters, V the applied voltage, F the Faraday constant, R the gas constant and T the absolute temperature; while η and λ describe transition rates for the transition between the open and the inactivated state (O2⇌I).

Defining as the probability of being in a specific state S at time t leads to the equation for the time evolution of the channels’ open probability , where the first two terms represent all transitions entering state O1 and the rightmost term all transitions leaving state O1. For sufficiently large numbers of the same channel, the quantities in this equation can be replaced by their macroscopic interpretation and the probability of being in a state S can be interpreted as the fraction of channels in S. The transition probabilities (transition rates times dt) can be described as rate constants, r, defining the number of channels changing from S to S in a given time period.[48] The total open probability ( number of open states), the ion channel number (Nc), the single channel conductance (g) and reversal potential of the ion (Eion) allow the calculation of the channels’ macroscopic current I(t,V) = Nc∙PO(t)∙g∙(V−Eion). The total membrane current results from the sum of the individual macroscopic currents considered and can be denoted as:

Since various electrophysiological phenotypes of the A549 cell could be predicted with high accuracy during the cell cycle, no additional leak current, summarizing all remaining channel activities, needed to be introduced.

In detail, HMMs were introduced in various studies [39,47,49-51] for the voltage-gated ion channels Kv1.3, Kv3.1, Kv3.4, Kv7.1 and CLC-2. The voltage-sensitive, but ligand-gated channel KCa1.1, pH and voltage-sensitive TASK-1 channels as well as the solely ligand-gated channel KCa3.1 were implemented according to Wang et al. [52], Limberg et al. [53] and Bailey et al. [54]. For CRACM1, a two state HMM (C⇌O) was defined and the corresponding rate constants for determining the open probability PO were derived from the measured, voltage-dependent closed and open lifetimes of this channel.[55] A constant current under consideration of the voltage-dependent conductivity was however implemented for the ligand-gated channels TRPC6 and TRPV3. The individual HMMs and corresponding rate constants of all ion channels are provided in S1 Fig and S1 Table.

The global cytosolic calcium concentration measured in A549 cells is about 64.7 ± 2.5 nM [56]. KCa3.1 (hIK) channels are activated at calcium concentrations greater than 200 nM and reach maximal activity at 1 μM.[57] Such high calcium levels required for activation indicate a close proximity of KCa3.1 to CRAC channels, where local calcium concentrations are much higher.[58,59] In the context of the calcium dependent gating of KCa3.1 and KCa1.1, the calcium inflow of CRACM1 channels is converted to a local calcium concentration by a formalism to depict the interaction between those channels in the model [39,60,61]: with where [Ca] is the intracellular calcium concentration, etrans a transfer coefficient, scaling the calcium influx within a certain intracellular subspace, in which the ion channels reside, to the calcium concentration, and ediff represents a calcium diffusion coefficient.[61]

Following Hou et al.[39], ediff was set to 3∙10−3 ms-1. The transfer coefficient etrans was estimated based on the cell volume, taking into account the mean measured cell capacitance Cmeasured = 33.4 pF and the specific A549 cell capacitance Cspecific = 2.45 μFcm-1 [62]. We assumed a cell compartment of 5% of the total cell volume (Volcompart = 2.336∙10−14 L), resulting in a transfer coefficient etrans of 21.9∙10−3 μMpA-1ms-1L-1.

CRAC channels are highly active at negative voltages, whereas their outward current, carried by Na+ at positive voltages, is negligible.[63] Changes in the calcium concentration provoked by CRAC channels occur slowly over time and thus take much longer than the applied test pulses. Hence, in order to reach the calcium concentration needed to trigger KCa3.1 and KCa1.1 activity, the steady state value after 10 seconds at a holding potential of -100 mV was taken as the input parameter for optimization and simulation. Nevertheless, the CRAC current evoked at negative voltages is considered in a time-dependent manner for the estimation of individual channel numbers and simulation of the initial and post pulse of performed patch clamp measurements. The time dependent evolution of the local calcium concentration is illustrated in S2 Fig.

Patch-clamp measurements of whole-cell membrane currents for model parametrization

A total of n = 16 of originally n = 50 A549 cell preparations, which met our internal lab standards (see Methods section), were considered for this study. A voltage-step protocol was applied consisting of an initial and re-pulse of -80 mV for 100 ms and a series of voltage pulses from -40 mV to +40 mV (increment 10 mV) of 800 ms duration. For model optimization, averaged whole-cell current curves were selected and used for fitting the whole-cell current simulated by the model. In addition, voltage-ramp measurements for determining the reversal potential and model verification were performed. The holding potential was set to -100 mV between all recording modalities. For detailed information on electrophysiological recordings, quality criteria and data pre-processing see “Methods” section.

Depending on the measured resting membrane potential and whole-cell current kinetics we were able to subdivide the cells into two groups, indicating the resting G0 and G1 phase of the A549 cell cycle. In detail, cells with a negative potential (Vrest = -40 mV to 0 mV, n = 11) showed a high instant activating current followed by a more slowly, time-dependent increase of the current (G0 phase). By contrast, cells with a highly depolarized, positive membrane potential (Vrest = 0 mV to +20 mV, n = 5) exhibited a comparable instantaneous current, but a lack of the time-dependent current, which resulted in lower steady currents (G1 phase) (see Fig 3A and 3B). The corresponding current-voltage curves were generated by selecting the steady state currents at 99% of the test pulse length. Averaged current-voltage curves with the corresponding standard deviations () for both cell groups are shown in Fig 3C and 3D. These variations can be explained by the expected heterogeneity in the morphology and membrane conductivity of the cell population.

Fig 3

Measured and simulated whole-cell current response.

Averaged whole-cell current curves measured at voltage-steps from -40 mV to +40 mV for (A) cells with negative resting membrane potential, G0 phase (n = 11) and (B) positive resting membrane potential, G1 phase (n = 5). (C) Averaged current-voltage curve of cells in G0 phase (mean ± SD): -40 mV: -0.1114 ± 0.0850 nA, -30 mV: -0.0849 ± 0.0671 nA, -20 mV: -0.0514 ± 0.0454 nA, -10 mV: -0.0059 ± 0.0292 nA, 0 mV: 0.0464 ± 0.0113 nA, +10 mV: 0.1137 ± 0.0286 nA, +20 mV: 0.1903 ± 0.0560 nA, +30 mV: 0.2770 ± 0.0920 nA, +40 mV: 0.3644 ± 0.1327 nA. (D) Averaged current-voltage curve of cells in G1 phase: -40 mV: -0.1140 ± 0.1721 nA, -30 mV: -0.0900 ± 0.1310 nA, -20 mV: -0.0649 ± 0.0960 nA, -10 mV: -0.0364 ± 0.0467 nA, 0 mV: -0.0065 ± 0.0088 nA, +10 mV: 0.0294 ± 0.0516 nA, +20 mV: 0.0744 ± 0.1150 nA, +30 mV: 0.1233 ± 0.1968 nA, +40 mV: 0.1685 ± 0.2794 nA. Comparison of simulated whole-cell current kinetics (blue lines) and experimental data (black lines) for the applied voltage-step protocol, and the corresponding current-voltage curves for cells in the G0 phase (E, G) and for cells in the G1 phase (F, H). RMSE values for simulation of activation curves are for the (E) G0 phase: RMSE = 0.0083 and (F) G1 phase: RMSE = 0.0074. RMSE values for current-voltage curves are for the (G) G0 phase: RMSE = 0.0031 and (H) G1 phase: RMSE = 0.0016. Averaged measured and simulated reversal potentials from derived current-voltage curves in the G0 phase are: Vrev_measured = -8.9 mV, Vrev_simulated = -9.3 mV (G) and in the G1 phase: Vrev_measured = 1.8 mV, Vrev_simulated = 0.95 mV (H). A comparison of the macroscopic ion currents in G0 phase and G1 phase can be found in S3–S7 Figs.

Fig 4 illustrates the measured resting membrane potentials, recorded reversal potentials from voltage-ramp measurements and reversal potentials from the generated current-voltage curves of the cells within the respective phases G0 and G1. Resting membrane potentials Vrest, reversal potentials Vrev determined by voltage-ramp measurements and derived from current-voltage curves differed statistically significantly between cells in G0 vs G1 phase. For detailed information on membrane potential measurement and statistics see “Methods” section and S2 Table.

Fig 4

Resting potential, ramp potential and reversal potential in G0 and G1 cell cycle phase.

Resting membrane potentials measured in current-clamp mode, reversal potentials from voltage-ramp measurements and derived reversal potentials from current-voltage curves in negative (G0 phase) and highly depolarized (G1 phase) A549 cells. Black lines in the boxplot indicate the median, white lines the mean value. Whiskers denote 10% and 90% percentiles. Single scores outside the 10% and 90% percentile are depicted as black circles. P-values demonstrate statistical significance between cells in G0 and G1 phase. Detailed measurement results are provided in S2 Table. Created with BioRender.com.

Model optimization and parametrization for simulation of the current-time characteristics based on experimental data

The HMM-based whole-cell model was implemented in the simulation environment MATLAB (MathWorks Inc.). Since HMMs model the open and closing probabilities of single ion channels stochastically, we instead used the macroscopic interpretation of the open probabilities PO to simplify model optimization in a deterministic manner and thereby reduce computational cost. Based on these open probabilities (PO), the macroscopic currents of all selected ion channel types were optimized by estimating the channel numbers () for the measured whole-cell currents using a particle swarm optimization (PSO) based approach:

This leads to an ordinary bounded integer least squares problem defined as [64]: where is the average measured total membrane current, is the matrix of the individual open probabilities and single channel conductivities , a vector containing the number of channels per channel type x and L and U are the lower and upper bounds of the individual channel numbers, respectively.

The model parameters are denoted in Table 1, showing the channel conductance (g), reversal potential of ions (Eion), bounds on the channel number for optimization and the estimated channel numbers for each single ion channel. The reversal potentials were calculated by the Nernst equation from the internal and external solutions used. Based on our experimental data and the available literature data, constraints for an adequate fitting of the measured whole-cell current curves were taken into account and set as follows: The minimal number of channels was set to 5 for Kv3.1, Kv3.4, Kv7.1, KCa1.1 and 10 for TASK-1, TRPC6, TRPV3 and CLC-2. Upper bounds were set to 150 for Kv3.1, 100 for Kv3.4 and CLC-2, 1350 for Kv7.1, 50 for KCa1.1 and 20 for TRPC6 and TRPV3. For Kv1.3 an amount of 20% and 80% [65] of the instant current were assumed as lower and upper bounds. The activity of KCa3.1 is strongly related to the cell cycle with which the current varies between 30% and 90% of the instant current depending on the membrane potential [66]. Cells with a hyperpolarized membrane potential exhibit a higher activity compared to cells with a depolarized membrane potential.[66] To account for the cell cycle dependent hIK activity in A549 cells, we set the lower bound for model optimization in the G0 phase to 50% of the instant current. For the G1 phase, 30% and 90% of the instant current were assumed as lower and upper bounds. In Leithner et al. [67], the TASK-1 current was determined to be about 80 pA (voltage level +40 mV) in hyperpolarized cells, estimated by 50 channels in the TASK-1 model. Thus, for TASK-1 the lower and upper bounds were set to 10 and 100 channels. The number of CRAC channels was set to 200 [68].

Table 1

A549 model parameters and estimated ion channel numbers from model optimization.

Ion channel	Single channel conductance	Eion	Constraints		Estimated channel number
			Lower	Upper	Negative cells G0—phase	Positive cells G1—phase
Kv1.3	15 pS [39,70]	-77.4 mV	20% of instant currenta: 15 (G0), 20 (G1)	80% of instant currenta: 59 (G0), 80 (G1)	22	20
Kv3.1	40 pS [49]		5	150	78	90
Kv3.4	14 pS [71]		5	100	5	54
Kv7.1	3.2 pS [72]		5	1350	1350	558
KCa1.1	250 pS [52]		5	50	40	15
KCa3.1	11 pS [39,57]		50% and 30% of instant currenta: 77 (G0), 63 (G1)	90% of instant currenta: 139 (G0), 188 (G1)	77	63
TASK-1	16 pS [73]		10	100	19	10
CRACM1	24 fS [68,74]	+95.6 mV	200		200	200
TRPC6	35 pS [75]		10	20	17	12
TRPV3	48 pS (negative voltages), 197 pS (+40 mV to +80 mV) [76,77]		10	20	12	13
CLC-2	2.8 pS [78]	-7.9 mV	10	100	13	11

a) At voltage-level +40 mV.

The measured whole-cell current curves (Fig 3A and 3B) obtained from voltage-step protocols at all voltage steps (-40 mV to +40 mV) were used for optimization. To optimize the number of ion channels instead of the least square error (see Eq 4), the relative root-mean square error (RRMSE) between the simulated (Imodel) and measured current (Idata), divided by the root mean square distance of Idata to a zero current trace, was minimized using particle swarm optimization [69]: where

Optimization by fitting the simulated curves to the experimental data was estimated using the averaged root mean square error values (RMSE) over all voltage steps:

Optimization results of the current time series data of cells assumed to be in the G0 and G1 cell cycle phase are shown in Fig 3E and 3F. The specified constraints (lower and upper bounds of the channel numbers) allowed a precise simulation of the current-time characteristics of the two phases, demonstrating an almost perfect fit (RMSEG0 = 0.0083, RMSEG1 = 0.0074) of the simulated whole-cell current kinetics (blue lines) to the experimental data (dashed black lines). In particular, the steady-state currents coincide well with the measured current amplitudes as shown in the derived current-voltage curves (Fig 3G and 3H). The source code of the A549 model for optimization of ion channel numbers and simulation of whole-cell currents is provided in S2 Text. Individual macroscopic currents of estimated channel numbers in G0 and G1 phase are illustrated in S3–S7 Figs.

When comparing the estimated ion channel numbers between the hyperpolarized and depolarized group the main changes are given in a decrease of ion channels generally responsible for hyperpolarization of the A549 membrane potential, including KCa1.1, KCa3.1 and TASK-1.[66,67,79,80] These changes coincide with the known alterations of ion channel activity and their influence on the membrane potential, and thus confirm the plausibility of the model optimization for both cycle phases. In addition, also the Kv7.1 current is lowered due to the lack of the time-dependent increase of the current, which equally provokes a slight depolarization of the cell membrane. By contrast, the number of Kv3.1 channels and the instantaneous activating current of Kv3.4 channels is rated higher by the model for depolarized than for the hyperpolarized cells, whereby the number of Kv1.3 is marginally decreased by 2 channels and TRPC6 by 5 channels, corresponding to a 34% decrease of the TRPC6 current.

Model verification based on experimental and synthetic data

For model verification voltage-ramp protocols from patch-clamp experiments were simulated and compared with the measured whole-cell currents for cells in G0 (n = 9) and G1 phase (n = 4). As shown in Fig 5A and 5B, the simulated ramp currents approximate well with the measured currents (RMSEG0 = 0.0314, RMSEG1 = 0.0211), confirming the accurate and reliable estimation of ion channel numbers by model optimization. Measured deactivation protocols were also simulated, but not considered further for model verification because of the limited quantity and low quality of experimental data (see S11 Fig).

Fig 5

Model verification by experimental and synthetic data.

Voltage-ramp currents measured at -100 mV to -40 mV (black lines) for (A) cells in G0 (n = 9) and (B) G1 (n = 4) phase, and corresponding model simulations (blue lines, RMSEG0 = 0.0314 and RMSEG1 = 0.0211). (C, D) Comparison of the whole-cell current of simulated sample cells (blue lines, n = 100) with averaged whole-cell current curves from activation protocols (black line, mean ± SD) in G0 and G1 phase.

As an additional verification step, we generated a synthetic data set of 1000 sample cells with the ion channel composition delivered from model optimization. In comparison, the macroscopic channel currents were simulated based on stochastic open probabilities of the individual Markov models (see Methods section) to test whether the expected open probability is a sufficient parameter for model parametrization and whether the model also represents the natural gating behavior in a large population of cells. The simulated whole-cell currents for all sample cells lay within the range of standard deviation of measurements, confirming the validity of this approximation for model optimization. Fig 5C and 5D exemplarily illustrates the simulated whole-cell current of 100 randomly selected sample cells at a voltage-level of +40 mV in G0 and G1 phase. Simulation results of all other voltage levels can be found in S8 and S9 Figs.

In-silico simulation of the A549 cell cycle

Optical control in performed patch-clamp experiments enabled us to exclude all cells during cell division (G2/mitosis phase). In particular, cells already indicating two nuclei and cells that were greatly enlarged in the antecedent G2 phase could not be patched because of their very weak and unstable cell membranes. Thus, patch-clamp experiments for model parametrization and verification comprise only those cells which were unambiguously in the G0, G1 or S cell cycle phase. According to the rhythmic oscillation of the membrane potential during cell cycle progression, showing a characteristic depolarization in G1 phase and hyperpolarization in G0 and S phase we were able to characterize patched cells being depolarized in G1 phase and hyperpolarized either in G0 or S phase.[19]

Since the optimized model fits well with the depolarization of the cell membrane by specific ion channel inactivation (c.f. [66,67,79,80]) demonstrating particularly a reduction of the TASK-1 and hIK current, we can distinguish hyperpolarized cells in the resting G0 phase from those in the further progressed S phase. Fig 6 shows all the channels involved and their changes in activity during the transition from G0 to G1 phase as predicted by the in-silico model.

Fig 6

Schematic illustration of the simulated ion channel activity during the A549 cell cycle.

G0/G1 transition: modulation of ion channel activity predicted by the model approach. G1/S and S/G2 transition: simulation of ion channel activation and inactivation, responsible for A549 cell cycle progression. Calculated membrane potentials for G0, G1, S and G2 phase are illustrated by yellow dots: Vm_G0 = -10.4 mV, Vm_G1 = -1.33 mV, Vm_S = -13.3 mV, Vm_G2/M = -5.49 mV. The symbols + and - indicate activation and inactivation of ion channels.

To estimate and predict membrane potential changes during cell cycle progression, the membrane potential Vm was calculated by setting Eq 7 to zero: where

Thereby, the potential V, for which the whole cell current Iwhole_cell equals zero in steady state condition, represents the membrane potential of the cell.

In addition, the membrane potential Vm was simulated by numerically solving the differential Eq 7 for Vm, starting at +5 mV until steady state was reached (t = 10 s, dt = 5.10−7). Both approaches led to almost the same membrane potential values and comply with the extracted reversal potentials Vrev from derived current-voltage curves and performed voltage-ramp measurements as shown in Fig 4. The time course of simulated membrane potentials in the different cell cycle phases is provided in S12 Fig. A comparison of measured and calculated potentials is shown in S3 Table.

Known actuators in G1–S transition of A549 cells are Kv1.3, Kv7.1, hIK and TRPV3 (Fig 6, G1 phase).[66,81-83] In accordance with the literature, we increased the current of these four channels, leading to a hyperpolarization to -13.3 mV, conceivable to trigger S phase initiation. In turn, the transition from S to G2/M phase again requires a depolarization of the membrane potential which is linked to an activation of TRPC6 channels in A549 cells.[84] A further activation by increasing the number of TRPC6 channels in the model results in a depolarization of the membrane potential to -5.49 mV, enabling the triggering of S-G2/M transition in the A549 cell cycle. In addition, we increased the CLC-2 current because of a well-known contribution of the CLC channel activation in S-G2/M phase initiation, although this does not lead to a noticeable increase of the membrane potential to –5.495 mV (Fig 6, G2/M phase).[19] The activation and inactivation of all ion channels involved and so far confirmed from literature and our own experimental data clearly demonstrates the characteristic oscillation of the membrane potential during cell cycle progression by simulations of the A549 in-silico model. The individual channel numbers and membrane potentials calculated for the different cell cycle phases are summarized in Table 2.

Table 2

Ion channel numbers and simulated ion channel blockages during different cell cycle phases in the A549 in-silico model.

	G0 phase			G1 phase				S phase			G2/M phase
A549 ion channels	optimized number of channels	simulated		optimized number of channels	simulated			Adjusted number of channels	simulated		Adjusted number of channels	simulated
		hIK block	TRPC6 block		Kv7.1 block	Kv1.3 block	TASK-1 block		hIK block	Kv1.3 block		TRPC6 block
Kv1.3	22	22	22	20	20	blocked	20	45a	45	blocked	45	45
Kv3.1	78	78	78	90	90	90	90	90	90	90	90	90
Kv3.4	5	5	5	54	54	54	54	54	54	54	54	54
Kv7.1	1350	1350	1350	558	blocked	558	558	650a	650	650	650	650
KCa1.1	40	40	40	15	15	15	15	15	15	15	15	15
KCa3.1	77	blocked	77	63	63	63	63	210a	blocked	210	210	210
TASK-1	19	19	19	10	10	10	blocked	10	10	10	10	10
CRAC	200	200	200	200	200	200	200	200	200	200	200	200
TRPC6	17	17	blocked	12	12	12	12	12	12	12	20a	blocked
TRPV3	12	12	12	13	13	13	13	20a	20	20	20	20
CLC-2	13	13	13	11	11	11	11	11	11	11	60a	60
Vm [mV]	-10.4	-3.96	-34.0	-1.33	+8.7	-1.15	+1.42	-13.3	+12.4	-13.0	-5.49	-31.1
ΔVm [mV]		+6.44	-23.6		+10.0	+0.18	+2.75		+25.7	+0.30		-25.6

a) Assumed changes of ion channel numbers correspond qualitatively to alterations of ion channel activity in the different cell cycle phases as described in the literature.

Model simulation and validation of altered ion channel activity of the A549 cell cycle based on literature data

Modulation of the membrane potential in different cell cycle phases, leading to promotion or interruption of cell cycle progression due to significant hyper- or depolarization, can be provoked by targeted activation and inactivation of ion channels.[19] Model simulation therefore provides an excellent tool for predicting the effects of ion channel activity on cell cycle progression. As a first validation step, we simulated 5 different, experimentally confirmed scenarios of ion channel activation and inactivation over the A549 cell cycle and their impact on the membrane potential and cell cycle progression.

A proven A549 cell cycle modulator is the TASK-1 channel. Its inhibition leads to a depolarization of the cell membrane and is related to reduced proliferation, mitosis and enhanced apoptosis.[67] Reduction or even a block of TASK-1 channel activity in G1 phase, as illustrated in Fig 7A, induces a depolarization of the membrane potential up to +1.42 mV, conceivable to arrest cells in G1 phase and thus prevent proliferation.

Fig 7

Schematic illustration of model simulations of ion channel activation and inactivation during cell cycle progression.

(A) Simulation of TASK-1 channel blockage in G1 phase. (B) Simulation of Kv1.3 channel blockage in a) G1 phase and b) during transition from G1 to S phase. (C) Simulation of Kv7.1 channel blockage starting at G1 phase. (D) Simulation of KCa3.1 channel blockage during G1-S transition. (E) Simulation of TRPC6 channel blockage in a) G0 phase and b) during transition from S to G2/M phase. The symbol + means activation, - means inactivation of ion channels.

Inhibition of Kv1.3 channels was shown to be associated with a depolarization of the cell membrane potential, accompanied by cell cycle arrest by impeding G1-S transition in A549 cells. Extending beyond this, recent studies have revealed reduced proliferation and suppressed tumor growth of A549 lung adenocarcinoma in vivo induced by Margatoxin (MgTX), a specific Kv1.3 channel blocker. The anti-proliferative effect of MgTX is related to proteins modulating the cell cycle, i.e. increased expression of the CDk1 inhibitor p21WAF1/Cip1 and a decrease of CDk4 and cyclin D3.[65,82] Simulating the inhibition of the Kv1.3 current in the G1 phase results in merely a small additional depolarization of the membrane potential from -1.33 mV to -1.15 mV (ΔVm = +0.18 mV) (Fig 7B, course a). Additionally, the blocking of these channels also results in a less negative membrane potential of -13.0 mV (c.f. unblocked state of -13.3 mV, ΔVm = +0.30 mV), while the KCa3.1, Kv7.1 and TRPV3 currents for G1-S transition remain unchanged (Fig 7B, course b). This lowered hyperpolarization is probably not sufficient to prevent G1-S transition itself, but the missing Kv1.3 current, respectively the changed membrane potential might affect the expression of cell cycle specific proteins, necessary for initiation of the G1-S transition.

Expression of Kv7.1 (KvLQT1) in A549 cells correlates with wound healing, motility and progression.[81] Blocking of the channels reduces wound healing rates by up to 31%, indicating a decrease of cell motility depending on the blocker concentration.[81] In addition, inhibition of Kv7.1 also significantly interferes with cell cycle progression by decreasing the proportion of cells in S and G2/M phases, resulting in a reduction of cell growth.[81] Simulation of a complete block of Kv7.1 channels in the G1 phase confirms these observations, leading to a significant depolarization up to +8.7 mV (Fig 7C), which is very likely to induce an overall cell cycle arrest.

Activation of hIK channels is known to promote G1/S transition in A549 cells, whereas inhibition increases the proportion of cells in G0/G1 phase. Model simulation of hIK channel inhibition in the G1 phase leads to a strong depolarization of Vm up to +12.4 mV (Fig 7D). This might suppress the transition from G1 to S phase and thus inhibit further cell cycle progression, resulting in a cell cycle arrest in G1 phase.[66,80,85] Similarly, the simulation of a knockdown or silencing of the channels by the model, this being equivalent to a total lack of hIK current starting in G0 phase, leads to a significant depolarization of the potential (Vm = -3.96 mV), possibly interrupting further progression. In addition, the depolarization in the G0 phase of ΔVm = +6.44 mV correlates with recent experimental data, confirming the estimated ion channel numbers for hIK by model optimization.[66]

TRPC6 channels are highly expressed during S-G2/M transition. Inhibition of these channels results in arrest of the cell cycle and decreased mitosis, invasion and proliferation.[84,86] Simulation of TRPC6 channel block in the S phase leads to a strong hyperpolarization of the estimated membrane potential (Vm = -31.1 mV), enabling suppression or bypassing of the required depolarization for S-G2/M transition and feasibly interrupting cell cycle and inhibition of mitosis as shown in Fig 7E, course b. By contrast, the knockdown of TRPC6 channels causes a further severe hyperpolarization of the cells in G0 phase (Vm = -34.0 mV), possibly to prevent the cells from entering into the cell cycle (Fig 7E, course a). The ion channel numbers and corresponding ion channel blockage of each simulation are given in Table 2.

Discussion

In this work we introduced for the first time an electrophysiological in-silico model of a cancer cell, in particular of the A549 human lung adenocarcinoma epithelial cell, taking into account the most relevant ion channels contributing to the membrane current and resting membrane potential of the cell. Single channel kinetics were modeled by a HMM approach and optimized by fitting the macroscopic currents of various ion channels functionally expressed in A549 cells, using whole-cell patch-clamp experiments. For model optimization, a linear optimization method (lsqlin Solver, Mathworks Inc.) and two non-linear, global optimization approaches (Simulated Annealing SA, simulannealbnd function Mathworks Inc. and Particle Swarm Optimization PSO, particleswarm function Mathworks Inc.) were used and the results compared. The PSO method showed, in comparison to the other algorithms, the most precise fit and highest stability and reproducibility over 1000 conducted simulation runs and was finally selected for model parametrization. For more information on the different optimization algorithms see S10 Fig.

The model takes into account the specific calcium and voltage dependencies of channel kinetics and enables the simulation and characterization of channel interaction, activation and inhibition, and prediction of membrane potential changes over parts of the cell cycle. The membrane potential of the cell is a fundamental physiological parameter, modulating various basic cellular functions, in particular proliferation of cells through rhythmic oscillation from hyperpolarized to depolarized states. These oscillations are caused by activation and inactivation of individual ion channels and trigger the transition between cell cycle phases [10,16,17,19].

The presented model facilitates a methodologically reliable and physiologically reasonable explanation of ion channel modulations and their impact on the whole-cell current and membrane potential as demonstrated in five cellular mechanisms published on the A549 cell line. The investigated mechanisms could be electrophysiologically explained and confirmed based on the simulated membrane potential changes. Specifically, inhibition of each of the channels Kv1.3, TASK1, KV7.1 and KCa3.1 induces a depolarization of the membrane potential and cell cycle arrest in G1 phase. Consistent with the literature, the cell cycle-specific block of each of these channels in the model in G1 phase provokes a depolarization of the membrane potential, plausibly representing the measured effect on the cell cycle. In comparison, inhibition of TRPC6 and activation of hIK channels results in a hyperpolarization of the cell membrane, impeding the S-G2/M phase transition and arresting the cells in the S phase. Simulation of both, TRPC6 channel block and increased hIK channel activity in S phase, led to a hyperpolarization of the membrane potential, also reliably reflecting the experimental findings from the literature. Beyond this, the absolute change in membrane potential when blocking the hIK channel activity in G0 phase is in accordance with experimental data, confirming the accuracy of optimization and validity of the model at this stage. The model therefore offers a first, valuable approach for investigating the effects of ion channel blockers or activators on the cell cycle and subsequently on cancer progression.

Nevertheless, a comprehensive experimental validation of this initial in-silico model of the A549 cell line is the next essential step. This includes a stepwise functional characterization of each single ion channel considered in the model by combined cell cycle analysis and patch-clamp experiments. Based on these extensive experimental investigations, the model can be further adapted and improved in order to finally provide a validated and powerful in-silico tool for research in cancer electrophysiology.

As in any electrophysiological cell model, the complexity of the biological system and limited experimental data require certain model assumptions and abstractions, which may restrict the accuracy and validity of such a model. In this approach, functionally expressed ion channels such as acid sensitive sodium channels (ASICs), amiloride sensitive sodium channels (ENaCs) and the CFTR chloride channels, regulated by cAMP, are not included in the current model. Also, Kv3.3, Kv9.3 and K2P9.1 channels were excluded because of a lack of available kinetic information. Nevertheless, the model contains all main functional ion current types and their specific kinetics present in A549 cells, comprising the large outward current sustained by potassium channels and, in comparison, the comparatively small inward current by chloride and calcium channels. Above all, however, the involvement of the main current types ensures a reliable and valid electrophysiological cell model, as demonstrated for a number of whole-cell ion-current models of excitable (cardiac or neural) cells at different levels of complexity in the past, allowing precise simulations of the action potential and its modulation, considering only a few summative ion currents [36-38,60].

As stated, current knowledge concerning ion channel expression and its function is limited for the modeled cell line. Indeed, the functional expression of more than 80% of ion channel genes found in A549 (95 of 113 channels) has not been proven yet. Thus, several ion channels were not considered in the current model state, limiting an accurate estimation of the remaining individual ion channel numbers by model optimization. The constraints for estimating the numbers of the respective ion channels in a realistic manner, however, were defined based on the experimental data from previous work and the literature data. Not only expression levels, but also the expression of a channel itself can basically be linked to distinct phases of the cell cycle, which represents a further uncertainty that cannot be assessed in the current state. For instance, a characteristic Cav3.1 current could only be found in 8% of A549 cells, which might be an indication of a cell cycle dependent expression of these channels in the plasma membrane.[87] Presently not considered ion channels as well as cell cycle specific expression levels can be attributed to deviations of the fitted whole-cell currents, in particular in the dynamics of the instant current at lower voltage levels (see Fig 3E and 3F), somewhat limiting the predictive power of the model in its current stage.

It is known that the kinetics of ion channels is strongly influenced by the biochemical environment and thus may differ between the different cell types.[88] Markov models available from other, non-excitable and transfected cell types, including HEK293, lymphocytes and oocytes, were used for the model approach which might affect the estimation of ion channel numbers in the model accordingly. In particular, the large number of Kv7.1 channels needed to simulate the time dependent increase of the whole-cell current in the G0 phase can be explained by alterations in the kinetics of the HMM model used, derived from Xenopus oocytes, assuming 70% of the channels in the inactivated state.[47] Similarly, the small deviations of simulated instant currents might be the result of slightly aberrant kinetics in the Markov models used. In addition, some of the parameters, such as single-channel conductivities, the calcium diffusion or transfer coefficient have not yet been quantified in this cell line and therefore needed to be adopted from other cell types respectively estimated for A549 cells. However, with regard to the calcium modeling it is important to note that the estimation of the local calcium concentration in this first model version is based on a simplified approach, assuming a certain steady state concentration which does not account for the complex local calcium dynamics.[61] Nevertheless, the local calcium distribution and time evolution in microdomains of the cell is highly significant for the activity of responding ion channels and, subsequently, affect the entire electric properties of the cell membrane.[89-92] Thus, to overcome this limitation and to further increase the model accuracy, a spatio-temporal modeling approach, addressing the heterogenous calcium profile and dynamics in the ER-PM junction provoked by local calcium release of CRAC channels from the edoplasmatic reticulum needs to be considered in future work of a more advanced model version.

The simulated and calculated membrane potentials Vm are based on the model parametrization of the measured activation curves, corresponding well to the measured reversal potentials Vrev from voltage-ramp protocols, but differ partly from the resting membrane potentials Vrest. The latter were directly measured after establishing the whole-cell configuration in current-clamp mode. In general, the measured membrane potentials (e.g. Vrest and Vrev) are strongly influenced by the experimental conditions. In whole-cell patch-clamp experiments the cytoplasm and the intracellular pipette solution are in an exchange process, which slightly alters the driving force over time and thus the membrane potential determined via the subsequently performed activation and ramp protocols.[93] Despite the convergence to zero, however, the reversal potentials derived from voltage-ramp measurements and the calculated membrane potentials demonstrate a significant difference between the hyperpolarized and depolarized cell groups in G0 vs. G1 phase. This difference is sufficient to confirm the possibility of membrane potential changes and their directions of change.

In summary, despite of these limitations, the introduced A549 model 1.0 allows a precise and reliable simulation of ion channel mediated alterations of the membrane potential in different cell cycle phases, serving as a first in-silico tool that supports the understanding of cancer electrophysiology in the human lung adenocarcinoma. Integration of continuously growing knowledge of ion channel expression and function in the A549 cell, together with further experimental investigations of ion channel expression, kinetics and function, and membrane potential changes over the A549 cell cycle will further increase the validity and predictive power of the model. We believe that this pioneering work may serve as a starting point for more advanced models in cancer electrophysiology taking into account additional biological, biochemical and electrophysiological information. Hodgkin-Huxley introduced the first mathematical whole-cell model that described the ionic mechanisms underlying the initiation and propagation of action potentials in an excitable nerve cell in 1952. Now almost 70 years later we have a first electrophysiological model of a cancer cell in hand, with the potential to usher in a new era of computational cancer electrophysiology.

Methods

A549 cell line

A549 human lung epithelial carcinoma cells, initiated from male lung carcinomatous tissue, were provided by the Center for Medical Research (ZMF) (Medical University Graz, Austria) obtained from the American Type Culture Collection (ATCC) and cultivated in Dulbecco’s Modified Eagle Medium (DMEM) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin. Cells were maintained at standard conditions in an incubator, 37°C, 5% CO2 in humidified atmosphere. The cells were split using Trypsin/EDTA every other day at confluence level of 50% to 70%.

Electrophysiological recordings and quality assessment

A549 cells (passage 6 to 8) were seeded on 6x6 mm coverslips in 6 well plates and cultured as described for 24h to 48h before electrophysiological recordings. Patch-clamp measurements were performed using the Axopatch-200B amplifier (Molecular Devices, USA) equipped with an inverted microscope (IM35, Zeiss, Germany) and Scientifica PatchStar micromanipulator. Patch-clamp pipettes were pulled from borosilicate glass capillaries (Assistant Mikro-Haematokrit capillary tubes, Hecht, Fisher Scientific) with a final resistance of 1.8 to 2.5 MΩ.

The external bath solution consisted of (in mM) 5.4 KCl, 137 NaCl, 5.6 d(+)-Glucose.H2O, 2.2 NaHCO3, 1.1 MgCl2.6H2O, 0.4 NaH2PO4.H2O, 10 HEPES/Na+, 1.8 CaCl2.2H2O, buffered with 2M NaOH to pH: 7.4. The intracellular solution contained (in mM) 4.3 ATP/K+, 110 KCl, 11 EGTA/Mg2+, buffered with KOH to pH: 7.4.

Resting membrane potentials were measured immediately after establishing the whole-cell configuration under current-clamp for 10 s and the entire trace averaged using Clampfit 10.3 software (Molecular Devices, USA). The resting membrane potential was low pass filtered at 50 Hz and traces digitized at 1 kHz using the Digidata 1322 A interface (Molecular Devices, USA) and Clampex 9.2 software (Molecular Devices, USA). Whole-cell currents were measured under voltage-clamp with a pulse protocol consisting of an initial- and re-pulse of -80 mV for 100 ms and 800 ms long test pulses from -40 mV to +40 mV (increment 10 mV). Voltage-ramp protocols for determining the reversal potential and model verification were performed starting from -100 mV to +60 mV in 20 s. The applied deactivation protocol consisted of an initial- and re-pulse of -40 mV for 150 ms, a depolarization pulse at 40 mV over 5 s for activation followed by seven 5 s long deactivation pulses from -40 mV to -100 mV (increment 10 mV). Currents were digitized with a sampling rate of 20 kHz and filtered at 2 kHz (Bessel). The holding potential was set to -100 mV between all recordings. Measurements were performed at room temperature between 22°C and 24°C.

To ensure a stable cell condition and reliable results, the experiments considered for data analysis were limited to a test time no longer than 20 minutes after removing cells from the incubator. The quality criteria for recordings are met showing a seal resistance at least greater than 1 GΩ, access resistance below 20 MΩ and a membrane capacitance greater than 20 pF. Cell capacitance was estimated according to a standard protocol from the patch-clamp system using the membrane test. In a data pre-processing step the capacitive peak of the current curves at the beginning of the test pulse of whole-cell recordings considered for data analysis and simulation were eliminated.

Computational modeling

The A549 model was implemented using the simulation environment MATLAB (R2019b, Mathworks Inc.). Macroscopic currents of single ion channels were modeled by hidden Markov models (HMMs) and the number of ion channels were optimized by fitting the macroscopic current to the measured whole-cell currents by a particle swarm optimization (PSO) algorithm (swarm size: 500, number of iterations: 5000, function tolerance: 1.10−6) using the Global Optimization Toolbox (Mathworks Inc.). In order to obtain a more precise optimization and, if possible, to ensure a global solution, we additionally used the fmincon solver as hybrid function after PSO optimization. The constraints for optimization were gradually restricted further and the solver was run again until the function value (fval) finally reached a constant value (fval = 1.5412.10−17) over all simulation runs, indicating a stable solution of the estimated ion channel numbers. The individual Markov models and corresponding rate constants of ion channels used in the A549 whole-cell current model are denoted in S1 Fig and S1 Table. In accordance with the applied voltage-step protocol in patch-clamp measurements, a corresponding pulse protocol was implemented for simulations. To consider the holding potential between the single measurements and to bring the channels into a defined initial state, we added an additional initial pulse of -100 mV for 100 ms in the pulse protocol for simulations (see S1A Fig). All hidden Markov models were simulated with a sampling frequency of 2000 kHz (step size dt = 5.10−7).

The matrix below exemplarily shows the state model of the Kv7.1 channel for simulating the expected open probability of the channels for model parametrization. The vector , representing the fraction of channels in each of the states S = {C1, C2, O1, O2, I}, was updated for each time step with the lsim function (Mathworks Inc.) for all voltage levels over the entire simulation protocol. Based on the expected open probabilities ( and ), the number of ion channels () was optimized by fitting the resulting individual macroscopic current to the measured whole-cell current.

The expected single channel currents for simulation and optimization of the macroscopic currents are illustrated in S1 Fig.

To model the stochastic opening and closing of the single ion channels in order to simulate the current of the sample cells, the hmmgenerate function (Mathworks Inc.) was used to generate a random sequence of states, depicting the single channel activity based on the Markov model. The transition probability matrix T for the Kv7.1 Markov model is defined as follows:

States corresponding to the channel being open (states 3 and 4) are summarized over all individual channels at each time step, constituting the number of channels used for calculating the macroscopic current according to Eq 3.

The source code for simulation of HMMs and model evaluation is provided in S2 Text.

Statistical analysis

Measurements of n = 16 of originally n = 50 individual A549 cells fulfilled the quality standards of our lab and were considered for data analysis. The cells were subdivided in two groups representing cells in the G0 phase (n = 11) and G1 phase (n = 5). Model optimization is based on the averaged measured whole-cell currents of the two cell groups. Mean values and standard deviations () of corresponding steady state currents at all voltage-levels are shown in Fig 3C and 3D. Normal distribution of measured resting potentials, ramp potentials and reversal potentials were tested using the Shapiro Wilk and the Kolmogorow-Smirnow tests. Statistical significance of measured resting potentials between cells in G0 phase and G1 phase was tested by using a two-tailed Student-t test. Statistical significance of reversal potentials from voltage-ramp measurements and current-voltage curves between the two groups was tested with the Mann-Whitney-U test. A p-value below 0.05 was considered significant (see Fig 4). All statistical analyses were performed using MATLAB (MathWorks Inc.) and results are summarized in S2 Table.

Literature review of ion channels in the A549 cell line.

(DOCX)

Click here for additional data file.

Source code A549 in-silico model.

(DOCX)

Click here for additional data file.

Parameters of hidden Markov models used in the A549 whole-cell current model.

Rate constants for transitions between the states of all hidden Markov models. For calcium dependent transitions the calculated steady state calcium concentration Ca_i = 4.68 μM is used.

(DOCX)

Click here for additional data file.

Patch-clamp results and statistical analysis.

Resting potentials of current-clamp measurements, reversal potentials from voltage-ramp measurements (ramp potential) and reversal potentials derived from current-voltage curves (reversal potential) of the individual cells with negative resting potential in G0 and positive resting potential in G1 phase.

(DOCX)

Click here for additional data file.

Comparison of measured and calculated membrane potentials in the different cell cycle phases.

(DOCX)

Click here for additional data file.

Kinetic schemes and expected single channel currents of implemented ion channels.

(A) Voltage-step protocol of patch-clamp measurements and below the corresponding protocol for model simulation. Kinetic scheme and expected single channel currents for voltage levels between -40 mV to +40 mV of the ion channels (B) Kv1.3, (C) Kv3.1, (D) TRPV3, (E) Kv3.4, (F) Kv7.1, (G) TRPC6, (H) TASK-1, (I) CRACM1, (J) CLC-2, (K) KCa1.1 and (L) KCa3.1.

(TIF)

Click here for additional data file.

Time evolution of the local calcium concentration.

Simulated changes of the local calcium concentration provoked by CRAC channels at holding potential of -100 mV over 10 s. Starting point c[Ca2+]i = 0.0647 μM, steady state c[Ca2+]i = 4.6847 μM, etrans = 21.8976∙10−3 μMpA-1ms-1L-1, ediff = 3∙10−3 ms-1.

(TIF)

Click here for additional data file.

Macroscopic currents of single ion channels in G0 and G1 phase.

Macroscopic currents of single ion channels, estimated by model optimization in respect to whole-cell current for A: G0 phase and B: G1 phase at +40 mV.

(TIF)

Click here for additional data file.

Comparison of macroscopic currents of voltage-gated potassium channels in G0 and G1 phase.

Comparison of macroscopic currents of voltage-gated potassium channels Kv1.3, Kv3.1, Kv3.4 and Kv7.1 in G0 and G1 phase at voltage levels from -40 mV to +40 mV.

(TIF)

Click here for additional data file.

Comparison of macroscopic currents of potassium channels in G0 and G1 phase.

Comparison of macroscopic currents of potassium channels TASK-1, KCa1.1 and KCa3.1 in G0 and G1 phase at voltage levels from -40 mV to +40 mV.

(TIF)

Click here for additional data file.

Comparison of macroscopic currents of calcium channels in G0 and G1 phase.

Comparison of macroscopic currents of calcium channels TRPC6, TRPV3 and CRAC in G0 and G1 phase at voltage levels from -40 mV to +40 mV.

(TIF)

Click here for additional data file.

Comparison of the macroscopic CLC-2 chloride current in G0 and G1 phase.

Comparison of the macroscopic CLC-2 chloride current in G0 and G1 phase at voltage levels from -40 mV to +40 mV.

(TIF)

Click here for additional data file.

Simulated whole-cell current of sample cells for G0 phase.

Simulated whole-cell current of 100 sample cells at (A) -40 mV, (B) -30 mV, (C) -20 mV, (D) -10 mV, (E) 0 mV, (F) +10 mV, (G) +20 mV, (H) +30 mV, (I) +40 mV for G0 phase. Black lines show the averaged measured whole-cell cell currents, background indicates the corresponding standard deviations at all voltage levels.

(TIF)

Click here for additional data file.

Simulated whole-cell current of sample cells for G1 phase.

Simulated whole-cell current of 100 sample cells at (A) -40 mV, (B) -30 mV, (C) -20 mV, (D) -10 mV, (E) 0 mV, (F) +10 mV, (G) +20 mV, (H) +30 mV, (I) +40 mV for G1 phase. Black curves represent the averaged whole-cell currents, grey background indicates the corresponding standard deviations at all voltage levels.

(TIF)

Click here for additional data file.

Fitting results of the different optimization algorithms applied.

Fitting results (n = 100 simulation runs) of (A, B) lsqlin, (C, D) SA and (E, F) PSO for averaged whole-cell currents of cells in G0 and G1 phase.

(TIF)

Click here for additional data file.

Simulated deactivation curves.

Simulation of the deactivation protocols for (A) G0 (n = 5, RMSEG0 = 0.0754) and (B) G1 (n = 3, RMSEG1 = 0.0673) phase. Deviations can be explained by the small sample size and apparent leakage currents for strongly negative deactivation pulses below -80 mV.

(TIF)

Click here for additional data file.

Simulated membrane potential in different cell cycle phases.

Simulated membrane potential (starting point at 5 mV over 10 s (dt = 5.10−7)) for (A) G0 phase Vm = -10.398 mV, (B) G1 phase Vm = -1.258 mV, (C) S phase Vm = -13.2 and mV (D) G2/M phase Vm = -5.263 mV.

(TIF)

Click here for additional data file.