Probing the Putative Role of K ATP Channels and Biological Variability in a Mathematical Model of Chondrocyte Electrophysiology

Abstract

Background:

Chondrocytes produce the extracellular cartilage matrix required for smooth joint mobility. As cartilage is not vascularized, and chondrocytes are not innervated by the nervous system, chondrocytes are therefore generally considered nonexcitable. However, chondrocytes do express a range of ion channels, ion pumps, and receptors involved in cell homeostasis and cartilage maintenance. Dysfunction in these ion channels and pumps has been linked to degenerative disorders such as arthritis.

Methods:

Because the electrophysiological properties of chondrocytes can be difficult to measure experimentally, mathematical modeling can instead be used to investigate the regulation of ionic currents. Such models can provide insight into the finely tuned parameters underlying fluctuations in membrane potential and cell behavior in healthy and pathological conditions. In this study, we introduce an open-source, intuitive, and extendable mathematical model of chondrocyte electrophysiology, implementing key proteins involved in regulating the membrane potential.

Results:

Because of the inherent biological variability of cells and their physiological ranges of ionic concentrations, we describe a population of models that provides a robust computational representation of the biological data. This permits parameter variability in a manner mimicking biological variation, and we present a selection of parameter sets that suitably represent experimental data.

Conclusion:

Our mathematical model can be used to efficiently investigate the ionic currents underlying chondrocyte behavior.

Introduction

An important feature of vertebrate joints is the presence of articular cartilage, the connective tissue covering the end of the bones, greatly reducing the friction of bone-against-bone and facilitating movement. Cartilage is largely composed of proteoglycans and collagens, forming a meshwork of extracellular matrix (ECM) proteins that can withstand large mechanical load while also remaining flexible compared with bone tissue.¹ As cartilage has poor regenerative properties, joint disorders such as osteoarthritis (OA) are becoming an increasing health burden as the global population ages.

OA is a progressive degradation of the cartilage in joints and is estimated to affect, for instance, 12.1% of adults over the age of 60 in the United States.² The progression of cartilage degradation, and what triggers it, remains poorly understood, and there are no treatments available other than pain relief or joint replacement surgery.³

Cartilage is largely made up of ECM. The cells responsible for the synthesis and turnover of this ECM are known as chondrocytes, embedded in small clusters in cartilaginous tissue. These cells are not innervated by the nervous system and are thus generally considered nonexcitable. In addition, as cartilage is not vascularized, chondrocytes rely on diffusion to exchange nutrients and waste material.⁴ Their extracellular milieu is therefore often slightly hypoxic, and chondrocytes have low metabolic turnover rates and poor regenerative abilities.⁵ This makes them prone to degenerative disorders, including OA, in which the cartilage generally is degraded at a faster rate than it is synthesized by chondrocytes.⁶

Although chondrocytes are considered nonexcitable, they do express a range of ion channels and ion pumps involved in the regulation of intracellular ionic concentrations and, in turn, the membrane potential and downstream intracellular processes. However, the exploration of electrophysiological properties of chondrocytes via traditional methods has presented challenges due to their small size (of <10 μm), representing a technical challenge for electrophysiological experiments and implying large uncertainties in measurement results.⁷

Mathematical modeling and computational simulation can provide an alternative, readily accessible methods of studying the membrane potential and ion dynamics of chondrocytes, and the construction of the model itself can be beneficial in understanding a target system's behavior. The mathematical model we present here can be seen as a reduced and abstracted representation of what is currently known about chondrocytes. In addition, such models can be useful to develop hypotheses, to identify mechanistic details or knowledge gaps, and to guide experiments.⁸ Mechanistic models built on biophysical knowledge such as the models presented here have an increasing interest for pharmacological research, as these have proved useful for qualitative understanding of the underlying physiology and pathophysiology.⁹

Moreover, these models can find direct application in drug development pipelines, see Strauss et al.¹⁰ Finally, computational models are also essential tools to study inter- and intravariability of cell physiology. Studying sources of potential variability can provide better understanding of biological processes and aid in making meaningful predictions. One intuitive approach to investigate sources of variability is to create and calibrate a population of models, which is used in the present study.¹¹

Maleckar et al.¹² earlier introduced a first-generation mathematical model for the membrane potential of the chondrocyte as based on an intensive literature review and additionally supported by unpublished experimental data. Both time-dependent and time-independent ion channel currents, as well as ion pump and ion exchanger currents, are represented in the model, with a detailed description of the underlying model development and equation derivation. This first-generation chondrocyte model focused mainly on K⁺ currents identified within the chondrocyte cell membrane. Furthermore, Maleckar et al. presented a model extension by including updated representation for the chondrocyte Na⁺/K⁺ pump.¹³

The previously published models were created using MATLAB,^12,13 a common commercial platform/software for model development and simulation. However, as a first-generation computational model of chondrocyte electrophysiology, we believe it beneficial to present our model on an open-source platform, so that more users can access it freely. Python is an open-source programming language that has been embraced by researchers across various branches of science, owing to its powerful libraries such as Numpy¹⁴ and Scipy¹⁵ in addition to interactive environments (IPython, Jupyter) that enable developers to provide interactive manuals of the software.

In this study, we present a reimplementation of the mathematical model for the chondrocyte's membrane potential first published in 2018 by Maleckar et al.¹² Our Python implementation is freely available at https://github.com/mmaleck/chondrocyte.

In this study, this model of chondrocyte electrophysiology is expanded by including the ATP-sensitive K⁺ (K_ATP) channel, a subgroup of inward-rectifying K⁺ channels that have been identified in human chondrocytes.¹⁶ This channel has been found to play an important protective role against hypoxia in the cardiovascular and nervous system,^17–19 particularly notable as oxygen tension in cartilage is lower than in vascularized tissue. K_ATP channel function is sensitive to intracellular Mg²⁺ concentrations, and channel activation is blocked by the binding of ATP.²⁰

Thus, when energy expenditure exceeds energy storage in the cell, intracellular ATP concentration decreases and the block of K_ATP channels is relieved,²¹ providing a direct link between cellular respiration and the membrane potential. The addition of K_ATP channel activity to the membrane dynamics of the model indicates that K_ATP channels may have a stabilizing effect on overall K⁺ currents in the presence of high Mg²⁺ concentrations.

To learn more about the model's behavior and parameter sensitivity, a population of models was created, allowing the investigation of ranges of parameter values and their effect on the steady state of the system. This investigation further clarified the important role of the Na⁺/K⁺ pump current in the modulation of the ionic balance of the model, as well as setting the resting membrane potential. Furthermore, the model can be expanded upon by the addition of ion channel and ion pump dynamics not described here. This readily available and expanded model of the chondrocyte thus provides insight into the dynamics of its membrane potential, and how this may be affected by variations in ion channel and ion pump function.

Methods

To investigate the chondrocyte membrane potential, we adopted a chondrocyte model and simulation protocols from Maleckar et al.^12,13 A schematic of the model is shown in Figure 1.

FIG. 1.

Schematic of the model and its implemented ion channels and ion pumps. INa−b, ICl−b, and IK₋b: background ionic activity for Na⁺, Cl⁻, and K⁺, respectively, INaCa: Na⁺/Ca²⁺ exchanger, INaH: Na⁺/H⁺ exchanger, INaK: Na⁺/K⁺ pump, ICa,ATP: ATP-dependent Ca²⁺ pump, IK−DR: delayed-rectifier K⁺ channel, IK−2p: two-pore K⁺ channel, IK−Ca: Ca²⁺ activated K⁺ channel, IK−ATP: ATP-sensitive K⁺ channel, ITRPV₄: Transient Receptor Potential V-4 (a mechanosensitive cation channel). Adapted from Maleckar et al.¹²

The given model is defined by a set of ordinary differential equations. Those systems are defined as $\begin{matrix} - d x ∕ d t \\ = f (x (t), u (t), k) \end{matrix}$ (1)

and describe the temporal evolution of states x(t). In the model used here, the states are the membrane potential (Vm) and the ion concentrations of Na⁺, K⁺, and Ca²⁺. Furthermore, the function f depends on currents u(f) for included ion pumps, channels and exchangers, and model parameters k. Systems described in the form as given in Eq. (1) need to be solved numerically.

The equation, which governs the chondrocyte transmembrane potential, follows the Hodgkin–Huxley formalism given in Eq. (2), where each cell component is represented as an electrical element, where cell membranes have a capacitance Cm and ion channels are treated as resistors. This gives a mathematical formalism for the membrane potential Vm, which includes all transport processes that are electrogenic and reads as follows: $C_{m} * \frac{d V_{m}}{d t} = \sum_{i o n s}^{} - I_{i o n}$ (2)

This work presents a model extension by adding the K_ATP current to the previously published chondrocyte model. Furthermore, a population of models¹¹ is used to investigate model behavior rigorously under a variety of different parameter distributions. The role of the maximum Na⁺/K⁺ pump current has particularly been studied in depth.

Mathematical modeling of the ATP-sensitive K⁺ current

We have added an ATP-sensitive K+ channel to the model from Maleckar et al.¹³ to investigate its potential functional role in regulating the membrane potential of chondrocytes.

The mathematical formulation of the K_ATP current IK−ATP reads: $I_{K} - A T P = σ g_{0} p_{0} f_{A T P} (V_{m} - E_{K})$ (3)

where σ is the channel density, g₀ is the unitary channel conductance, p₀ is the maximum open channel probability, and f_ATP is the fraction activated channels. E_K denotes the equilibrium potential for K⁺ in the given circumstances. In the previous implementation of IK−ATP, a constant value was used for g₀. However, the unitary conductance can be expressed as²² $g_{0} = γ_{0} f_{M} f_{N} f_{T}$ (4)

where γ₀ is the unitary conductance in the absence of intracellular Na⁺ and Mg²⁺ and depends on [K⁺]₀: $γ 0 = 35.375 {(\frac{{[K^{+}]}_{0}}{5.4})}^{0.24}$ (5)

The term f_M in Eq. (4) represents the inward rectification generated by intracellular Mg²⁺ ions and can be expressed by means of a Hill equation: $f_{M} = \frac{1}{1 + \frac{{[M g^{2 +}]}_{i}}{K_{h, M g}}}$ (6)

where K_h,_Mg is the half-maximum saturation constant that depends on membrane potential and on [K⁺]₀:

${K_{h M g}}^{0} ({[K +]}_{0}) = 0 . 65 ∕ s q r t ({[K +]}_{0} + 5)$ (8)

Here, the value of the electrical distance δ_Mg was set to 0.32, F is the Faraday constant, R is the gas constant, and T is the absolute temperature. In a similar manner, the term f_N represents the inward rectification caused by intracellular Na⁺ ions and is expressed as $f_{N} = \frac{1}{1 + {(\frac{{[N a^{+}]}_{i}}{K_{h, N a}})}^{2}}$ (9)

K_{h, N a} = K_{h, N a}^{0} e x p (- \frac{δ_{N a} F}{R T} V_{m})

(10)

where δ_Na = 0.35 and K_h⁰_,Na = 25.9.

Population of models

As introduced in the prior section, in the context of our work, we use the term population of models for set model simulations with randomly varied parameter sets for ionic current conductances.¹¹ Generally, a population of modes is useful to investigate parameter sensitivity,^23,24 sources of variability,²⁵ as well as emergent model behavior dependency for a wide range of parameters.

In this study, we searched for parameters that had a significant effect on the simulation results and the physiological relevance of the model. Table 1 offers an overview of the selected set of parameters, the values in the original implementation, the parameter regime used to create a population of models, and a brief parameter description.

Table 1.

Ensemble of Parameters Selected for Creating a Population of Models

Parameter	Initial value	Selected range	Parameter description
Na_i	25 (mM)	17.5–32.5 (mM)	Na⁺ initial concentration
Ki	180 (mM)	126–234 (mM)	K⁺ initial concentration
Cai	1 (μM)	0.75–3.25 (μM)	Ca²⁺ initial concentration
INaK,scale	1.625	1.375–2.1125	Scaling factor: max Na⁺/K⁺ pump current
gNab,bar	0.1 (pS)	0.07–0.13 (pS)	Background Na⁺ leakage conductance

Parameter ranges originate from parameter sampling from log normal distribution. The parameter distribution was not fitted to experimental data and therefore underlies the assumptions that the parameters' distributions are skewed and the parameter values are positive. The distribution underlies the reasoning that numerous works have shown that a log normal distribution is often a useful assumption to describe the random variation in biological samples.²⁶ The population of models presented in the results is created out of 100 simulation runs. Each member of the population has its characteristic set of parameters, and all the parameter sets were sample from log normal distributions with the mean being the parameter value and variance σ = 0.15.

Results

The following results section is divided into three parts. We first provide, in brief, illustrative model validation for our new implementation that is used for the following sections. In the second part, we present simulations focused on the K_ATP currents and the effect of different Mg²⁺ concentrations. Finally, the parameter and model behavior is studied, based on the population of models approach.

Validation

We performed the validation against a previous publication, Maleckar et al.¹³ Figure 2 shows an exemplar: the results of temperature-dependent contribution of the Na⁺/K⁺ pump electrogenic current to the chondrocyte resting membrane potential, as implemented both in MATLAB and Python, replicating Figure 2 from Maleckar et al.¹³ Our model includes many other currents (as expounded in the original publication), all compared side-to-side with one another to ensure that the new Python implementation was faithful to the original implementation.

FIG. 2.

Temperature-dependent contribution of the Na⁺/K⁺ pump electrogenic current to the chondrocyte resting membrane potential at (A) 23°C and (B) 37°C. Unbroken lines indicate our new implementation in Python, while dotted lines indicate the published MATLAB implementation.

ATP-sensitive K⁺ current

To investigate the effect of IK−ATP on K⁺ dynamics in the chondrocyte, numerical simulations were performed, as shown in Figure 3. Accurate measurements of intracellular Mg²⁺ concentrations in chondrocytes are, to our knowledge, unavailable, hence, we tested a range of initial Mg²⁺ values. The Mg²⁺ concentrations used for the simulation were (a) 0.1 mM, (b) 1.0 mM, and (c) 10 mM, while three different K⁺ concentrations, 5, 30, and 70 mM, were used for each Mg²⁺ concentration as indicated in the Figure 3. To further clarify the effect of different intracellular Mg²⁺ concentration to the overall chondrocyte matrix, we introduce Figure 4, which (a) shows the steady state voltage dependence of IK−ATP while (b) illustrates the time-dependent changes of chondrocyte membrane potential under different intracellular Mg²⁺ concentrations.

FIG. 3.

Sensitivity of IK−ATP against varying extracellular K⁺ concentration based on different intracellular Mg²⁺ concentrations. Intracellular Mg²⁺ concentrations used are (A) 0.1 mM, (B) 1.0 mM, and (C) 10 mM, respectively, while three concentrations of extracellular K⁺ concentrations (5, 30, and 70 mM) are used for all panels.

FIG. 4.

Intracellular Mg²⁺ concentration-dependent contribution of the ATP-sensitive K⁺ current to the chondrocyte resting membrane potential. In (A), we illustrate steady-state voltage dependence and the Mg²⁺ dependence of the IK−ATP current, while (B) shows the hyperpolarization of chondrocyte membrane potential with different intracellular Mg²⁺ concentration. The extracellular K⁺ concentration used for (A, B) is fixed at 7 mM.

Populations of models

An example population of 100 models with randomly varied sets of parameters (Table 1) is illustrated in Figure 5. The steady state solutions for all species are reached within the simulation time. Notably, the steady states for both Na⁺ and K⁺ reach similar concentration levels regardless of the set of parameters, whereas Ca²⁺ concentrations at steady state solution are strongly affected by the set of parameters. Thus, the behavior of the membrane potential follows the shape of the Ca²⁺ time curve, as influenced by the composition of the parameter sets.

FIG. 5.

Population of 100 models with randomly drown scaling factor for the maximum of the sodium-potassium (Na⁺/K⁺) pump current and background sodium (Na⁺) leakage conductance. Each simulation trajectory results from a parameter set with unique parameters. Steady-state values for the resting membrane potential, V, (A) and the intra-cellular calcium concentration, [Ca²⁺]_I (C), are affected by the variation in these parameters, as is the time point at which intracellular potassium and sodium concentrations, [K⁺]_I (B) and [Na⁺]_I (D), reach steady states. Each simulation trajectory results from a randomly drawn parameter set; a sampling profile for the example parameter Ca_i,0 is shown in inset. A simulation time of 2500 s ensures that all species reach steady state.

For all 100 parameter sets, the Na⁺ steady state solution is close to zero. However, the time point when reaching the steady state is parameter-dependent (Figs. 5B, D and 6A, B). To investigate the individual parameter effects, simulations with randomly drawn parameters for each individual parameter represented in Table 1 were performed. Larger parameter ranges were also investigated; in Figure 5 INaK,scale ranges from 1.375 to 2.1125, whereas for Figure 6, INaK,scale values between 0.1 and 7.0 were tested.

FIG. 6.

The scaling of the maximum sodium-potassium (Na⁺/K⁺) pump current, INaK,scale, affects the timing of the intracellular sodium concentration Na+, [Na+]i, steady state as well as the steady state of the chondrocyte resting membrane potential, V. The simulation trajectories for V and [Na+]i in the low (A) and the middle (B) parameter regime of INaK,scale are displayed here. Simulations for INaK,scale > 5.0 cannot be presented as negative intracellular Na+ concentrations cause the simulation to fail. A and B display simulation curves for 0.1 < INaK,scale < 0.6. The steady state for [Na+]i ranges from 10 mM up to 90 mM. B and D show simulation trajectories for 0.6 < INaK,scale < 5.0. In all cases, Na+ depletion in [Na+]i can be observed. The timing is parameter-dependent.

The simulation results indicate that the scaling factors of the maximum Na⁺/K⁺ pump current (INaK,scale) has a major effect on the membrane potential, the final concentrations for Na⁺ and K⁺, as well as the duration required for the system to reach steady-state. A wide range of INaK,scale values resulted in a Na⁺ depletion, which is unlikely to occur in a physiologically viable chondrocyte preparation. A stepwise scan through INaK,scale (0.1 < INaK,scale < 7.0) reveals three parameter regimes for the scaling of this current. For 0.1 < INaK,scale < 0.6, here called the low parameter regime, Na⁺ concentrations range from 10 to 90 mM (Fig. 6A). If 0.6 < INaK,scale < 5.0, a nonphysiological state of Na⁺ depletion is eventually reached, as intracellular sodium moves toward zero over time.

The value of INaK,scale also affects the time course of this depleted state (Fig. 6B): the smaller the scaling factor and smaller the current, the later the system reaches its steady state (the middle parameter regime). For INaK,scale > 5.0, steady-state Na⁺ concentrations become negative, which is biologically and numerically infeasible and causes the simulation to fail.

A population of models approach was again used to further probe the model dependence on parameter variation for low and middle parameter regimes, this time varying initial conditions (Fig. 7). Figure 7 shows strong model dependence on initial conditions for evolving variables (panels A–D). However, while there are a variety of stable steady states at different parameter sets, where for example, intracellular sodium is not depleted (panel B), all these occur for the low parameter regime of INaK,scale (blue traces, Fig. 7, all panels).

FIG. 7.

Study of interplay between the initial conditions and I_NaK,scale on the simulation outcome for membrane potential (A), intracellular sodium concentration (B), the intracellular potassium concentration (C), and the intracellular calcium concentration (D). Fifty parameter sets with randomly drawn initial conditions for evolving concentrations and membrane potential (Na_i_0, Na_o, K_i_0, K_o_0, Ca_i_0, Ca_o, and V_0; initial values for the intracellular sodium concentration, the extracellular sodium concentration, the intracellular potassium concentration, the extracellular potassium concentration, the intracellular calcium concentration, the extracellular calcium concentration, and the resting membrane potential, respectively) were created. Simultaneously, 10 different values for INaK,scale through the low to middle parameter regimes (0.1–4.9, 0.5 step) were used, for a total of 10 × 50 simulations. The colors reveal this INaK,scale gradation, with dark blue as the lowest value (0.1), and deep red as the highest value (4.1). The same colors always represent the same INaK,scale value, but different overall parameter sets.

Thus, the low parameter regime of $I_{N a K}$ conductance (0.1 < INaK,scale < 0.6), corresponds to, roughly, an ${\bar{I}}_{N a K}$ current of approximately 0.9–5.35 pA/pF and likely represents a physiologically relevant model regime incorporating a variety of stable steady states for evolving ionic concentrations and the membrane potential. To further investigate the overall dependence of the chondrocyte's resting membrane potential on K⁺ currents, we varied the conductance parameters for these at a variety of INaK,scale values in the low parameter regime (Fig. 8). While there are a variety of stable steady states for the resting membrane potential of the chondrocyte in this example, even as INaK,scale increases (moving right from panel A to C), all the values for each unique parameter set in the current model set the chondrocytes resting membrane potential somewhere between approximately −50 and −70 mV.

FIG. 8.

The resting membrane potential evolution within the low parameter regime for INaK,scale, as this increases, as the background potassium conductance (g_K_b_bar) is varied via the population of models approach. Simulations (100 parameters sets, varying I K,b: g_K_b_bar, the conductance governing the background potassium current, over a distribution as described in Methods) were run with differing values for INaK,scale in the low parameter regime range (0.1–0.6) until intracellular sodium reached steady state. Each color in panels (A–C) above refers to a unique parameter set, and panels (A–C) show the chondrocyte resting membrane potential reaching steady state over time for increasing values of I_NaK,scale (0.1, 0.3, and 0.5, respectively).

Discussion

In this study, a model of chondrocyte electrophysiology was expanded to include the ATP-sensitive K⁺ (KATP) channel and an extended ionic “population of models” approach was used to understand quantitatively the steady-state behavior of the chondrocyte membrane potential.

ATP-dependent potassium channel and the chondrocyte model

One possible interpretation of Figure 3 is that higher Mg²⁺ concentrations help to stabilize IK−ATP currents with respect to varying extracellular K+ concentrations. The lower Mg²⁺ concentration additionally contributes to the hyperpolarization of the membrane potential (Fig. 4). This stabilizing effect is interesting, as experimental data suggest a protective role of high Mg²⁺ on joint health. Although the exact mechanisms behind this protection remain elusive, some hypotheses point to a role of Mg²⁺ in reducing chondrocyte apoptosis and facilitating chondrocyte proliferation.^27–30

Previous work has also shown that the sulfonylurea glibenclamide stimulates growth of human chondrocytes by insulin-like growth factor 1 (IGF-1)-dependent mechanisms.^31,32 As IGF-1 levels are inversely related to inflammatory markers and oxidative stress, and promoted by nutrients, including Mg²⁺,³³ it bears speculation that Mg²⁺ regulatory effects on chondroprotection are multifold. One caveat of the present study is that physiological Mg²⁺ concentrations in chondrocytes remain unknown, and thus we have tested our model with a range of concentrations. Reliable measurements of intracellular Mg²⁺ concentrations will be necessary to assess the physiological effects of Mg²⁺ on K⁺ currents, and the potentially protective role of Mg²⁺ and the ATP-dependent K⁺ channel in cartilage maintenance.

Insights from the population of models approach

To get a complete understanding of the three parameter regimes for the scaling factor of the maximum Na⁺/K⁺ pump, current further investigation is needed. Notably, we have found that in the low parameter regime for this scaling factor, a plethora of stable steady states for the chondrocyte model is reached, setting the emergent chondrocyte resting membrane potential somewhere between approximately −50 and −70 mV. However, this investigation has not taken the variety of potential isoforms of the three subunits underpinning the Na-K-ATPase and its generated current, as done in our previous work.¹³ In addition, it would be interesting to further investigate the link between Ca²⁺ concentration and the membrane potential in a subsequent study.

Limitations of the study and possible model improvement

The current-voltage relationship for IK-ATP for differential extracellular potassium concentrations as shown in Figure 2 follows that of Zhou et al,³⁴ the parent formulation used in this implementation. As in the original, the reversal potential (Nernst potential) for potassium was fixed at −94 mV in this case. However, in the context of widely fluctuating or increased extracellular potassium concentration, which is a distinct possibility in the environs of the chondrocyte, it is likely that the Nernst potential will influence the driving force for the current, and should be calculated explicitly, to be further considered in subsequent model improvements.

Furthermore, the IK-ATP model used here is neither specific to a chondrocyte nor indeed to a pure isoform, but rather to models of the mammalian cardiac sarcolemma. Furthermore, while the present work has appropriately considered how IK-ATP is modulated by changes in intracellular magnesium, other important ionic processes, also known to be significantly changed when intracellular magnesium is altered, have not yet been considered in this model iteration, including the sodium potassium pump, the 2-pore potassium channel, and the calcium pump, among others. Future model development will address this gap.

In addition, an important present limitation is the relative importance of the plasma membrane versus the mitochondrial IK-ATP channel, important in neuro- and cardioprotection. Here, we have only accounted for the surface membrane channel and have not yet implemented mitochondrial regulation. Even if K-ATP channels were unchanged at the level of the plasma membrane, these might still exert effects in the intracellular domain.

Clearly, the lack of available biological data presents a challenge for the validation of models such as the one presented here. While investigating ranges offer insights into plausible parameterization, measurement data are also needed to circumscribe these ranges directly. Due to the high parameter uncertainty of the model, care must be taken in interpreting results to make useful predictions and to suggest hypotheses for further testing. A more thorough investigation into the range of appropriate intracellular and synovial ionic concentrations remains to be performed.

Concluding Remarks

Here, we present a model of chondrocyte electrophysiology, with particular attention to parameter sensitivity and the role of K⁺ currents in setting the steady-state membrane potential. Our preliminary results support a putative role for Mg²⁺ in cartilage maintenance as recorded in clinical studies, by potentially stabilizing K⁺ currents through the ATP-dependent K⁺ channel, although more research is required to bolster this speculation.

Furthermore, we show how a population of models can be used to examine fluctuations in a range of parameters, and their interactions as the model reaches steady-state membrane potential. Such random variations in parameter values offer a more robust representation of reality: naturally occurring fluctuations and uncertainty are inherent in biological, including electrophysiological, data. Finally, the chondrocyte model has been reimplemented from MATLAB into Python to increase its accessibility and is now available as an open-source repository on GitHub, with demo scripts to aid interested parties to get started immediately.

Footnotes

Acknowledgments

The authors acknowledge the members of our research teams and collaborators for their support and encouragement.

Author Disclosure Statement

No competing financial interests exist.

Funding Information

The authors gratefully acknowledge funding for this work from the Norwegian Ministry of Education and Research via the Simula-UiO-UCSD Research PhD Training initiative (SUURPh Programme) as well as the Research Council of Norway, which funds PROCardio, a Norwegian Centre for Research-Based Innovation (SFF).

References

Woo

, Buckwalter

. Injury and repair of the musculoskeletal soft tissues. J Orthop Res, 1988; 6:907–931.

Dillon

, Rasch

, Gu

, et al. Prevalence of knee osteoarthritis in the united states: Arthritis data from the third national health and nutrition examination survey 1991–1994. J Rheumatol, 2006; 11:2271–2279.

Grässel

, Muschter

Recent advances in the treatment of osteoarthritis [version 1; peer review: 3 approved. F1000Research, 325(9(F1000 Faculty Rev)):1–17, 2020.

Archer

, Francis-West

. The chondrocyte. Int J Biochem Cell Biol, 2003; 4:401–404.

Rajpurohit

, Koch

, Tao

, et al. Adaptation of chondrocytes to low oxygen tension: Relationship between hypoxia and cellular metabolism. J Cell Physiol, 1996; 2:424–432.

Akkiraju

, Nohe

. Role of chondrocytes in cartilage formation, progression of osteoarthritis and cartilage regeneration. J Dev Biol, 2015; 4:177–192.

Wilson

, Clark

, Banderali

, et al. Measurement of the membrane potential in small cells using patch clamp methods. Channels (Austin), 2011; 5:530–537.

Wolkenhauer

Why model?. Front Physiol, 2014; 5:21.

van der Graaf

, Benson

, Peletier

. Topics in mathematical pharmacology. J Dyn Differ Equ, 2016; 28:1337–1356.

10.

Strauss

, Wu

, Li

, et al. Translational models and tools to reduce clinical trials and improve regulatory decision making for QTc and proarrhythmia risk (ich e14/s7b updates). Clin Pharmacol Ther, 2021; 109:319–333.

11.

Britton

, Bueno-Orovio

, Van Ammel

, et al. Experimentally calibrated population of models predicts and explains intersubject variability in cardiac cellular electrophysiology. Proc Natl Acad Sci, 2013; 110:E2098–E2105.

12.

Maleckar

, Clark

, Votta

, et al. The resting potential and K+ currents in Primary Human articular chondrocytes. Front Physiol, 2018; 9:1–21.

13.

Maleckar

, Martín-Vasallo

, Giles

, et al. Physiological effects of the electrogenic current generated by the Na+/K+ pump in mammalian articular chondrocytes. Bioelectricity, 2020; 2:258–268.

14.

Harris

, Millman

, van der Walt

, et al. Array programming with NumPy. Nature, 2020; 585:357–362.

15.

Virtanen

, Gommers

, Oliphant

, et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat Methods, 2020; 17:261–272.

16.

Mobasheri

, Gent

, Nash

, et al. Evidence for functional ATP-sensitive (K_ATP) potassium channels in human and equine articular chondrocytes. Osteoarthr Cart, 2007; 15:1–8.

17.

Yamada

, Ji

, Yuan

, et al. Protective role of ATP-sensitive potassium channels in hypoxia-induced generalized seizure. Science, 2001; 29:1543–156.

18.

Crawford

, Jovanovic´

, Budas

, et al. Chronic mild hypoxia protects heart-derived H9c2 cells against acute hypoxia/reoxygenation by regulating expression of the SUR2A subunit of the ATP-sensitive K+ channel. Memb Transp Struct Funct Biogen, 2003; 278:31444–31455.

19.

Turzo

, Vaith

, Lasitschka

, et al. Role of ATP-sensitive potassium channels on hypoxic pulmonary vasoconstriction in endotoxemia. Respir Res, 2018; 19:29.

20.

Trapp

, Tucker

, Ashcroft

. Activation and inhibition of K-ATP currents by guanine nucleotides is mediated by different channel subunits. Proc Natl Acad Sci USA, 1997; 16:8872–8877.

21.

Noma

ATP-regulated K+ channels in cardiac muscle. Nature, 1983; 305:147–148.

22.

Ferrero

, Sáiz

, Ferrero

, et al. Simulation of action potentials from metabolically impaired cardiac myocytes: Role of ATP-sensitive K+ current. Circ Res, 1996; 79:208–221.

23.

Marder

, Taylor

. Multiple models to capture the variability in biological neurons and networks. Nat Neurosci, 2011; 14:133–138.

24.

Vagos

, Arevalo

, de Oliveira

, et al. A computational framework for testing arrhythmia marker sensitivities to model parameters in functionally calibrated populations of atrial cells. Chaos, 2017; 27:093941.

25.

Mould

, Upton

. Basic concepts in population modeling, simulation, and model-based drug development—part 2: Introduction to pharmacokinetic modeling methods. CPT Pharmacomet Syst Pharmacol, 2013; 2:1–14.

26.

Limpert

, Stahel

, Abbt

. Log-normal distributions across the sciences: Keys and clues: On the charms of statistics, and how mechanical models resembling gambling machines offer a link to a handy way to characterize log-normal distributions, which can provide deeper insight into variability and probability—Normal or log-normal: That is the question. BioScience, 2001; 51:341–352.

27.

Feyerabend

, Witte

, Kammal

, et al. Unphysiologically high magnesium concentrations support chondrocyte proliferation and redifferentiation. Tissue Eng, 2006; 12:3545–3556.

28.

, Ma

, Pang

, et al. Synthesis of chondroitin sulfate magnesium for osteoarthritis treatment. Carbohydr Polym, 2019; 212:387–394.

29.

Yao

, Xu

, Zheng

, et al. Intraarticular injection of magnesium chloride attenuates osteoarthritis progression in rats. Osteoarthr Cartil, 2019; 27:1811–1821.

30.

Kuang

, Chiou

, Lo

, et al. Magnesium in joint health and osteoarthritis. Nutr Res, 2021; 90:24–35.

31.

Heinze

, Vetter

, Holl

, et al. Glibenclamide stimulates growth of human chondrocytes by IGF I dependent mechanisms. Exp Clin Endocrinol Diabetes, 1995; 103:260–265.

32.

Heinze

, Vetter

, Fussgänger

, et al. The sulfonylurea glibenclamide directly stimulates growth of chondrocytes. Pediatr Res, 1988; 23:112.

33.

Maggio

, De Vita

, Lauretani

, et al. IGF-1, the cross road of the nutritional, inflammatory and hormonal pathways to frailty. Nutrients, 2013; 5:4184–4205.

34.

Zhou

, Cortassa

, Wei

, et al. Modeling cardiac action potential shortening driven by oxidative stress-induced mitochondrial oscillations in guinea pig cardiomyocytes. Biophys J, 2009; 97:1843–1852.