Action potentials in vitro: theory and experiment

Action potential generation underlies some of the most consequential dynamical systems on Earth, from brains to hearts. It is therefore interesting to develop synthetic cell-free systems, based on the same molecular mechanisms, which may allow for the exploration of parameter regions and phenomena not attainable, or not apparent, in the live cell. We previously constructed such a synthetic system, based on biological components, which fires action potentials. We call it “Artificial Axon”. The system is minimal in that it relies on a single ion channel species for its dynamics. Here we characterize the Artificial Axon as a dynamical system in time, using a simplified Hodgkin-Huxley model adapted to our experimental context. We construct a phase diagram in parameter space identifying regions corresponding to different temporal behavior, such as Action Potential (AP) trains, single shot APs, or damped oscillations. The main new result is the finding that our system with a single ion channel species, with inactivation, is dynamically equivalent to the system of two channel species without inactivation (the Morris-Lecar system), which exists in nature. We discuss the transitions and bifurcations occurring crossing phase boundaries in the phase diagram, and obtain criteria for the channels’ properties necessary to obtain the desired dynamical behavior. In the second part of the paper we present new experimental results obtained with a system of two AAs connected by excitatory and/or inhibitory electronic “synapses”. We discuss the feasibility of constructing an autonomous oscillator with this system.

1 Introduction

The physics of excitable media is largely concerned with the patterns in space and time created by nonlinear excitations. In this context, electrophysiological processes, and action potentials specifically, hold a unique place as some of the most consequential dynamical systems on earth. Action potentials (APs) play a vital role in biological computation, as sequences of APs encode information in a variety of ways [1].

We have recently introduced a minimal synthetic system, the “Artificial Axon” (AA) [2–4], which is capable of generating APs in time. The experimental system is based on a traditional suspended lipid bilayer (“black lipid membrane”) with embedded voltage gated potassium ion channels (KvAP). An ionic gradient maintained across the membrane provides the free energy source to elicit action potentials. Non-traditionally, the system is held in the off-equilibrium, excitable state by a modified voltage clamp (“Current Limited Voltage Clamp”: CLVC) which allows for voltage dynamics. The system is minimal in that it is built with one ion channel species only, yet it can support APs. The key is the dynamics of the channel, which includes inactivation. Previously we reported on the dynamics of firing APs in the AA. We showed that the threshold behavior of the system is the same as in real neurons [5], namely it is governed by a saddle node bifurcation [4].

Here we first discuss numerically the phase diagram for a single AA, identifying regions in parameter space which give rise to trains of APs, damped oscillations, or single shot APs. We propose a simplified version of the Hodgkin-Huxley model [6] for our system, based on the kinetics of our voltage gated ion channel KvAP [7]. The main new insight is that a system with one ion channel species, with inactivation, is dynamically equivalent to the (biological) system with two channel species, without inactivation (the Morris-Lecar system [8]). This latter system has been studied extensively in real biological contexts. In their seminal paper on the barnacle giant muscle fiber, Morris and Lecar display experimentally the variety of behavior the barnacle axon can produce, including AP trains and damped oscillations, and confirm that it is mainly due to two voltage dependent “conductances”, corresponding to $K^{+}$ and $C{a}^{++}$. They then show, through a 3D dynamical systems model “of the Hodgkin-Huxley form”, that indeed a system of two voltage gated, non inactivating conductances, with no further elements, can generate most of the behavior observed in the experiments [8]. Comparisons between the AA and the Morris-Lecar system will punctuate our presentaton, as the dynamic behavior of the two systems is similar. We identify the bifurcations separating different regions in the phase diagram of the AA; in particular, we point out a transition which may not have been described before in elctrophysiology. In the second part of the paper we present experimental measurements of the inactivation and recovery rates of the channels in the AA system. These help to qualitatively place the present experimental system with respect to the phase diagram obtained from the model, and understand the requirements on channel dynamics in order to obtain autonomous oscillations in the AA. Finally, we demonstrate a system of two AAs connected by electronic “synapses”, as a prototype for future network developments.

2 Theory

2.1 The artificial axon system

In the Artificial Axon, the phospholipid bilayer acts in essence as the dielectric of a parallel plates capacitance, sandwiched between two conducting media which are the electrolytes on either side. This capacitance is charged by two kinds of ionic currents: the current through the ion channels embedded in the membrane, and the current sourced by the clamp electrodes. The charge carriers for the former are K⁺ ions, for the latter Cl⁻, Ag⁺, and all other ions in solution. The voltage dynamics $V(t)$ is governed by the following Equation 1 [4]:

$\begin{array}{c}\hfill \frac{dV}{dt}=\frac{N_0\chi }{C}{p}_o\left(t\right)\left[V_N-V\left(t\right)\right]+\frac{{\chi }_c}{C}\left[V_c-V\left(t\right)\right]\end{array}(1)$

where $C$ is the membrane capacitance, $N_0$ the number of ion channels, $\chi$ the open channel conductance, $p_o(t)$ the fraction of channels in the open (conducting) state (so $N_0\chi p_o(t)$ is the total channel conductance); $V_N$ is the Nernst potential corresponding to the bulk concentrations of $K^{+}$ ions on the two sides of the membrane, ${\chi }_c$ is the CLVC conductance, and $V_c$ the CLVC command voltage (which is a control parameter in the experiments). The first term on the RHS of Equation 1 is the channel current (divided by $C$), proportional to the driving force $[V_N-V(t)]$, since $V_N$ is the equilibrium potential and $V(t)$ the present potential. The second term is the clamp current, proportional to the driving force $[V_c-V(t)]$. This second term is exactly equivalent to the presence of a second, reversed ionic gradient (of sodium ions, say) with Nernst potential equal to $V_c$ and corresponding leak conductance ${\chi }_c$ [2, 4]. Equations of the basic form Equation 1 underlie many models of nerve excitability [9].

The probability that channels are open, $p_o(t)$, is determined by the voltage dependent channel dynamics. The minimal model for the KvAP channel has 3 states: open, closed, and inactive (Figure 1). In order to minimize the dimensions of parameter space, we connect them with 4 rate constants as shown in the figure. The unidirectional arrows are, strictly speaking, unphysical, but since here we are concerned with the macroscopic dynamics of the AA rather than the microscopic channel dynamics, they represent a permissible approximation. In fact, the detailed channel dynamics is more complex than shown in Figure 1, with more states and corresponding transition rates [7]. We assume that these complications do not change the qualitative features of the phase diagram of the system, nor the nature of the bifurcations which occur in the dynamics. The rate equations which determine $p_o(t)$ in Equation 1 are then:

$\begin{array}{c}\hfill \left\{\begin{array}{c}\frac{d{p}_i}{dt}=p_o\left(t\right)k_i-p_i\left(t\right)k_r\hfill \\ \frac{d{p}_o}{dt}=\left[1-p_o\left(t\right)-p_i\left(t\right)\right]k_o\left(V\right)-p_o\left(t\right)\left[k_c\left(V\right)+k_i\right]\hfill \end{array}\right.\end{array}(2)$

where $p_i(t)$ is the probability that the channels are in the inactive state $(p_o+p_i+p_c=1)$, and the rates are as in Figure 1. Since the channel is voltage gated, $k_o$ and $k_c$ are voltage dependent, which couples Equation 2 to Equation 1. We assume a standard Arrhenius dependence [7, 10]:

$k_o\left(V\right)=\kappa e^{\alpha \left(V-V_0\right)},k_c\left(V\right)=\kappa e^{-\alpha \left(V-V_0\right)}(3)$

this symmetric form being chosen once again to minimize the number of parameters [4, 8]. In the experimental system, the inactivation and recovery rates, $k_i$ and $k_r$, are similarly voltage dependent; however this dependance does not qualitatively alter the dynamical phases and bifurcations of the system. To simplify the discussion, we take these rates as constant throughout Section 2. We refer to the model (1), (2) with voltage independent $k_i$, $k_r$ as the “voltage independent model”.

Figure 1

Figure 1. Simplified model for the KvAP channel dynamics. The three states are open (O), closed (C), and inactive (I). Adapted from [4], licensed under Creative Commons Attribution License (CC-BY 4.0), doi:10.1088/2399-6528/ac43d0.

In order to understand the requirements on channel dynamics to obtain various temporal patterns, in Section 2.2 we construct a diagram in parameter space of the dynamical behavior of the system (1), (2), representing a single AA. We refer to it as a “phase diagram” for short. The difficulty is that even in this minimal representation, the parameter space is still high dimensional. To make progress, we identify the most relevant parameters in relation to the experiments: the clamp voltage $V_c$ and clamp conductance ${\chi }_c$, which are the actual control parameter in the experiments, and the effective rates of Figure 1, which define the suitability of the channel for obtaining interesting dynamical behavior, such as autonomous firing. In general (and specifically for the KvAP), the rate of opening and closing are much faster than those of inactivation and recovery, $k_o,k_c\gg k_i,k_r$, and so the most relevant parameters are then the clamp voltage, the clamp conductance [11], and the rates of inactivation and recovery.

2.2 Phase diagram

Figure 2A shows one representation of the phase diagram (in dimensional variables) for the system (1), (2) representing the AA, namely a cut through parameter space for constant ${\chi }_c=5\times 10^{-10}{\mathrm{\Omega }}^{-1}=500$pS and $V_c=-50$mV, with $k_i$ and $k_r$ voltage independent. The diagram is generated by simluating AAs with the given parameters and recording the resulting time traces.

Figure 2

Figure 2. (A) Phase diagram of the dynamic behavior obtained from the model (1), (2), with voltage independent inactivation and recovery rates $k_i$ and $k_r$. (Figure 1). The figure is a cut through a higher dimensional parameter space, at constant ${\chi }_c=500$pS and $V_c=-50$mV. Region I corresponds to AP trains. Region II exhibits smaller amplitude “oscillations”, distinct from APs. Region III corresponds to damped oscillations, and Region IV to single shot APs. (B) A “heat map” representing the same region in parameter space as (A), with the color indicating the frequency of oscillation (in mHz). The central region with high frequency oscillations indicate region II and III, while the darker regions on the left and right are region I and IV, respectively.

We identify four regions of distinct behavior. In region I the dynamical system (1), (2) produces action potential trains (i.e. limit cycles). The firing rate increases for increasing $k_r$, while the width of the AP decreases for increasing $k_i$. Figure 3 shows a corresponding time trace (blue). In region II the system exhibits “oscillations”, distinct from Region I in that the amplitude is smaller, and the rate higher. The transition from region I to region II can be sharp, depending on the control parameters ${\chi }_c$ and $V_c$. In Figure 3 we plot time traces on the two sides of the transition, for comparison. We explore this transition in further detail in the next section.

Figure 3

Figure 3. Two representative time traces of the voltage independent model, illustrating the sharp increase in frequency as one crosses from Region I to Region II in the phase diagram of Figure 2A. The blue trace has $k_r=0.18$${\mathrm{s}}^{-1}$ while the red trace has $k_r=0.19$s⁻¹, with all other parameters identical ($k_i=10.4$s⁻¹). Firing is elicited by stepping $V_c$ from −200 mV to −50 mV at $t=1$s. The purple and orange traces show the probability that channels are open, $p_o$, for the blue and red traces respectively, with scale on the second y-axis.

Region III corresponds to damped oscillations, the damping increasing for increasing $k_r$; representative time trace are shown in Figure 4 (red) and Figure 5 (blue). In Figure 4 we display the behavior of the time traces as we move across the transition $II\to III$. Finally, in region IV the system fires only once, after which the voltage remains constant at a relatively high value, corresponding to the limit cycle collapsing to a stable fixed point different from the resting potential. The red trace in Figure 5 displays this behavior.

Figure 4

Figure 4. Two time traces of the voltage independent model, displaying the transition from region II to region III in the phase diagram of Figure 2A. The blue trace ($k_r=0.2$s⁻¹) corresponds to region II, the red trace ($k_r=0.22$s⁻¹) to region III. Parameters other than $k_r$ are identical to Figure 3.

Figure 5

Figure 5. Two time traces of the voltage independent model, displaying damped oscillations and single shot AP. The blue trace ($k_r=0.22$s⁻¹) corresponds to region III of the phase diagram of Figure 2A, the red trace ($k_r=0.5$s⁻¹) to region IV. The purple and orange traces display the corresponding probability that the channels are inactive, $p_i$. Parameters other than $k_r$ are identical to Figure 3.

The lines separating the different regions in Figure 2A were obtained by qualitative assessment of time traces long enough to determine, for example, whether oscillations are damped or not. A more quantitative measure is displayed in Figure 2B, where the color indicates the firing rate. The yellow-orange region corresponds to regions II and III, while the blue region on the upper left is region I and the deeper blue part on the lower right is region IV. To summarize, for the case that the inactivation and recovery rates ($k_i$ and $k_r$) are voltage independent, we find four possible dynamic behaviors under constant current input: AP trains, oscillations, damped oscillations, and single shot firing. Remarkably, this phenomenology of the AA (one channel species with inactivation) is the same as for the Morris-Lecar system (two channel species without inactivation) [8, 12, 13].

To further characterize the nature of the transition between regions I and II of the phase diagram of Figure 2A, we look at the behavior of the firing rate across the transition. Figure 6A shows the firing rate vs. $k_r$ for various values of the clamp voltage $V_c$, with fixed $k_i=10.4$s⁻¹, and ${\chi }_c=500$pS. For $V_c=-54$mV we see a sharp transition at $k_r\approx 0.21$s⁻¹ where the firing rate increases steeply. This transition smoothes out as the clamp voltage is raised, with the transition point moving to lower values of $k_r$. Figure 6B displays the same transition for 3 different values of $k_i$, keeping $V_c=-50$mV fixed. We see that the main effect of varying $k_i$ is to shift the transition point.

Figure 6

Figure 6. (A) Firing rate as a function of $k_r$ for several values of $V_c$ (legend) with $k_i=10.4$s⁻¹. For each curve, the lower firing rate portion (to the left of the “corner”) corresponds to Region I of the phase diagram, while the higher firing rate portion corresponds to Region II. In this example, the transition is sharp for $V_c=-54$mV. As $V_c$ is raised, the sharpness of the transition decreases and the transition region moves to smaller values of $k_r$. For each curve, the last four to five points at the higher frequencies lie in Region III (damped oscillations). (B) Firing rate as a function of $k_r$, for different values of $k_i$ (legend) and fixed $V_c=-50$mV. This plot displays the same transition as in (A). For both plots, ${\chi }_c=500$pS.

The phase diagram of Figure 2A is a slice through a higher dimensional parameter space with $V_c$ and ${\chi }_c$ held fixed. Changing these parameters (within a reasonable range) shifts the boundary lines in the $k_i$, $k_r$ plane but does not fundamentally alter the possible behaviors (Supp. Mat. Supplementary Figures S2, S3). Taking a different slice in parameter space (e.g. ${\chi }_c$ vs. $V_c$) also yields the same regions of behavior (Supp. Mat. Supplementary Figures S4, S5).

2.3 2D system and bifurcations

Displaying phase space trajectories often gives better insight into the nature of a bifurcation, however the dynamical system (1), (2) is 3D, which makes phase space representations cumbersome. We can make progress by noting that, in the regime $k_i,k_r\ll k_o,k_c$, the 3D system can be reduced to 2D. We will further assume that the reduction does not affect the nature of the bifurcations [8]. We may also cast the problem in dimensionless variables: Equation 1 suggests the choice of $\tau =C/(N_0\chi )$ as the time scale and $V_N$ as the voltage scale; then the dimensionless control parameters in (1) are the clamp voltage $\widetilde{V_c}=V_c/V_N$ and clamp conductance $\widetilde{{\chi }_c}={\chi }_c/(N_0\chi )$, while the dimensionless rates are $\widetilde{k_i}=\tau k_i=k_{i}C/(N_0\chi )$ and similarly for ${\tilde{k}}_r$.

To reduce the system to 2D, we consider that interconversion between open and closed states of the channels is faster than all other time scales in the system; we may then assume that the balance $C\rightleftharpoons O$ has at all times its “equilibrium” value at the given voltage. We introduce a new coordinate $p_a(t)$, the probability that channels are “active” (i.e. not inactive: open or closed). Since the interconversion $C\rightleftharpoons O$ is “fast”, we may make the substitution $p_o(t)\to p_a(t)P(V)$ in (1), where $P(V)$ is the equilibrium opening probability in the absence of inactivation, which is a function of voltage only. In terms of rates, $P(V)=1/(1+k_c/k_o)$. The function $P(V)$ is measured in the experiments [14]; it is indeed well approximated by a Fermi-Dirac distribution (see (Equation 3)):

$P\left(V\right)=\frac{1}{k_c\left(V\right)/k_o\left(V\right)+1}=\frac{1}{e^{-\frac{2\alpha }{T}\left(V-V_0\right)}+1}(4)$

Writing the total channel conductance $N_0\chi$ as ${\chi }_0$ (so $\tau =C/{\chi }_0$), Equation 1 can be made dimensionless by dividing all voltages by $V_N$ and scaling by $\tau$. With the condition $p_a(t)+p_i(t)=1$, the system (1), (2) can now be written in terms of the coordinates ($V(t)$, $p_a(t)$) as:

$\begin{array}{c}\hfill \frac{dV}{dt}=p_a\left(t\right)P\left[V\left(t\right)\right]\left[1-V\left(t\right)\right]+{\chi }_c\left[V_c-V\left(t\right)\right]\end{array}(5)$

$\begin{array}{c}\hfill \frac{d{p}_a}{dt}=k_r-k_r\left[1+\frac{k_i}{k_r}P\left[V\left(t\right)\right]\right]p_a\left(t\right)\end{array}(6)$

where we have introduced dimensionless variables as mentioned above, dropping the tilde, so $V/V_N\to V$, $t/\tau \to t$, ${\chi }_c/{\chi }_0\to {\chi }_c$, etc. Numerically, using “standard values” for our experimental parameters [4], $C=300$pF, $\chi =200$pS, $N_0=100$, the time scale is $C/(N_0\chi )=1.5\times 10^{-2}$s, so for example the (dimensional) rate $k_r=0.2$s^-1 corresponds to the dimensionless rate $(C{k}_r)/(N_0\chi )=\tau k_r=3\times 10^{-3}$. If we regard the parameters ($\alpha$, $V_0$) which define the open probability function $P(V)$ as fixed, the 2D dynamical system (5), (6) depends only on the four control parameters ${\chi }_c$, $V_c$, $k_r$, and $(k_i/k_r)$.

In reducing the dynamical system from 3D to 2D, we used the general idea that “fast variables” can be considered to attain their “equilibrium value” for time scales over which the “slow variables” vary. This same procedure is used for example in [15], as well as by Morris and Lecar [8]. The assumption is that the reduction does not change the qualitative behavior of the system. While it is not easy to give a rigorous justification, below we show numerically that the 2D system indeed displays dynamic behavior corresponding to all four regions in the diagram of Figure 2A. Following [15], we show in Figure 7 time traces for the 3D and the 2D system, in Region I. As expected, there is qualitative similarity, but quantitative differences, in the amplitude and frequency of the oscillations, for the two models.

Figure 7

Figure 7. A visual comparison of AP trains in the 3D system (Equations 2, 3, blue) to AP trains in the 2D system (Equations 5, 6, red). The blue trace is the same as in Figure 3, while the red trace is generated by taking the same initial conditions as the blue trace ($V_c$ stepped from −200 mV to −50 mV at $t=0.2$ s) and converting to dimensionless units as described in the text, with $k_i=0.15$, $k_r=9.2\times 10^{-3}$. Both traces are in Region I of their respective phase diagrams, in the vicinity of Region II.

Referring to the phase diagram of Figure 2A as a guide, let us explore the nature of the transitions as we move along a horizontal line in that phase diagram (i.e. we vary $k_r$ at fixed $k_i$, other parameters fixed). For the fixed values $V_c=-1.7$, $k_i=0.15$, ${\chi }_c=0.05$ we find the behavior of region I (AP trains) in the interval $4.03\times 10^{-3}\lesssim k_r\lesssim 9.20\times 10^{-3}$. Figure 8A shows a phase space trajectory for $k_r=6.0\times 10^{-3}$ (blue line), which is a limit cycle. The corresponding AP train is shown in Figure 8B. In Figure 8A we also plot the nullclines $d{p}_a/dt=0$ (red line) and $dV/dt=0$ (orange line); they intersect at the unstable fixed point inside the limit cycle.

Figure 8

Figure 8. (A) Phase space trajectory (blue) in the $V$, $p_a$ plane for the 2D dynamical system (5), (6), displaying the limit cycle corresponding to an AP train, with $k_r=6.0\times 10^{-3}$, $k_i=0.15$, $V_c=-1.7$, and ${\chi }_c=0.05$, starting from the initial condition $(V=-1,p_a=1)$. Also shown are the nullclines $d{p}_a/dt=0$ (red) and $dV/dt=0$ (orange); they intersect at an (unstable) fixed point. (B) Time trace of the AP train corresponding to the limit cycle shown in (A); dimensionless $t$ and $V$ (see text).

For $k_r\ge 9.21\times 10^{-3}$ the system exhibits damped oscillations, corresponding to Region III in the phase diagram of Figure 2A. The transition between AP trains and damped oscillations is sharp. Figure 9 shows the same quantities as Figure 8, for $k_r=9.20\times 10^{-3}$ (just inside the region of AP trains). Figure 10 shows these plots for $k_r=9.21\times 10^{-3}$ (just outside the region of AP trains). Here the phase space trajectory spirals into the (stable) fixed point, corresponding to damped oscillations of $V(t)$. Note that with the above parameter values (specifically, the relatively small $k_i$) we did not cross Region II. Rather, we may say that the transitions $I\to II$ and $II\to III$ have “merged”, in the sense that the transition from AP trains to damped oscillations is also accompanied by a steep increase in frequency (Figures 9B, 10B). The phenomenology just described is that of a subcritical Hopf bifurcation [16]. The linear stability analysis of the fixed point close to the bifurcation shows that the eigenvalues of the stability matrix form a complex conjugate pair and cross the imaginary axis from right to left as the fixed point changes from unstable to stable. The eigenvalues $\lambda$ (real and imaginary part) are shown in Figure 11 for different values of $k_r$ near the bifurcation. With the parameter values of this section, the bifurcation point is at $k_r\approx 9.21\times 10^{-3}$.

Figure 9

Figure 9. (A) Phase space trajectory displaying the limit cycle for $k_r=9.20\times 10^{-3}$ (other parameters are same as in Figure 8), just inside Region (I) (B) Time trace of the AP train corresponding to the trajectory shown in (A).

Figure 10

Figure 10. (A) Phase space trajectory for $k_r=9.21\times 10^{-3}$ (other parameters are same as in Figure 8), just outside Region I. There is no longer a stable limit cycle and the trajectory spirals into the stable fixed point. (B) Time trace of the damped oscillations corresponding to the trajectory shown in (A).

Figure 11

Figure 11. Eigenvalues $\lambda$ of the stability matrix at the fixed point, calculated numerically for the dynamical system (5), (6), for different values of $k_r$. From left to right, the points correspond to: $k_r=(10.0,9.30,9.20,9.10,9.0,8.0)\times 10^{-3}$. Other parameters are as in Figure 8.

For $k_r\gtrsim 12\times 10^{-3}$ (not plotted) the system is so heavily damped that there are no oscillations; depending on initial conditions $V(t)$ either approaches the fixed point value from one side or fires once and then approaches the fixed point. This is the behavior of region IV of Figure 2A.

In general, the behavior of the system is independent of the choice of initial conditions. Since there is only one fixed point, for any starting point $(V(0),p_a(0))$ the system either spirals into the fixed point, if it’s stable, or moves to the limit cycle. However, there are exceptions to this near the Hopf bifurcation. Figure 12 shows such a case. The system starts at $(V,p_a)=(-0.5,0.39)$, which is close to the fixed point. The trajectory travels outwards and makes several loops before stabilizing at the larger stable limit cycle. This behavior arises from the presence of an unstable limit cycle which is in between the stable fixed point and stable limit cycle. For $k_r$ values far from the Hopf bifurcation, this phenomenon is not seen because the fixed point itself is unstable.

Figure 12

Figure 12. (A) Phase space trajectory for $k_r=9.18\times 10^{-3}$ and initial condition $(V=-0.5,p_a=0.39)$, showing the presence of an unstable limit cycle inside the outer stable cycle. Other parameters are as in Figure 8. (B) Time trace of the trajectory shown in (A).

Moving up in the phase diagram, i.e. increasing $k_i$, we recover Region II in the 2D system as well. This is displayed in Figure 13, which is analogous to Figure 6, but obtained for the 2D system. In Figure 13A we plot the firing rate vs. $k_r$ with $k_i$ fixed at $k_i=0.35$ and different values of $V_c$. It is evident that exactly the same phenomenology occurs as for the 3D system: at a critical value of $V_c$ the transition is sharp ($V_c=-1.73$ with these parameter values), and it smoothens out as $V_c$ is raised, with the transition shifting to smaller values of $k_r$. Figure 13B displays the same transition for $V_c=-1.7$ fixed and different values of $k_i$ (legend).

Figure 13

Figure 13. (A) Firing rate for the reduced 2D model as a function of $k_r$, varying $V_c$ (legend) with $k_i=0.35$ and $\chi =0.05$ held fixed. The phenomenology is the same as for the 3D system (Figure 6), in particular, there is a sharp transition for a critical value of $V_c$. (B) The same transition displayed for fixed $V_c=-1.7$, ${\chi }_c=0.05$ and different values of $k_i$ (legend).

In Figures 14, 15 we show representative phase space trajectories and time traces across the Region I $\to$ Region II transition, for the 2D system. Different from the Hopf bifurcation corresponding to the transition II $\to$ III, the fixed point inside the limit cycle remains unstable on both sides of the transition. This is confirmed by an analysis of the eigenvalues of the stability matrix at the fixed point, which remain on the right side of the imaginary axis in both cases.

Figure 14

Figure 14. Phase space trajectory (A) and corresponding time trace (B) for the 2D system just prior to the I $\to$ II transition. $k_r=17.03\times 10^{-3}$, $k_i=0.35$, $V_c=-1.73$, ${\chi }_c=0.05$.

Figure 15

Figure 15. Phase space trajectory (A) and corresponding time trace (B) for the 2D system just after the I $\to$ II transition. $k_r=17.05\times 10^{-3}$, $k_i=0.35$, $V_c=-1.73$, ${\chi }_c=0.05$

In addition to the subcritical Hopf bifurcation, the system also contains a few other bifurcations which are similar to those of the Morris-Lecar system [12, 13]. Using $k_i$ as the bifurcation parameter, we briefly describe some of the bifurcations encountered as $k_i$ is increased. Starting at small $k_i$, there is one stable fixed point at the intersection of the $V$ and $p_a$ nullclines which is globally stable. As $k_i$ increases, the V nullcline moves to the right, and a saddle-node bifurcation will occur when the two nullclines intersect at a second point, which splits into two additional fixed points post-bifurcation. As $k_i$ is further increased, one of the newly created fixed points will annihilate with the original stable fixed point in another saddle-node bifurcation. This causes limit cycles to arise (Region I and Region II), as the only remaining fixed point is unstable. Finally, as $k_i$ is increased further, the remaining fixed point becomes stable through the Hopf bifurcation described above and all trajectories spiral into it (Region III), until eventually no oscillations occur (Region IV).

The transitions described above may of course be explored following different trajectories in parameter space. Since for the 2D system (5), (6) the control parameters which have to do with channel rates are $k_r$ and $k_i/k_r$, a natural trajectory is to keep the latter fixed. The overall picture remains the same: as an example, we show in Figure 16 the firing rate vs. $k_r$ for fixed $k_i/k_r$; the different curves correspond to the clamp values in the range $-1.72<V_c<-1.52$, in increments of 0.02.

Figure 16

Figure 16. The transition from region I to region II in the 2D system, explored along trajectories with fixed $k_i/k_r=25$ (rather than fixed $k_i$ as in Figure 13C). The different curves correspond to $V_c$ in increments of 0.02, starting at $V_c=-1.52$ (violet) and ending at $V_c=-1.72$ (red).

2.4 Analogy to the magnetization transition

The plots of Figures 6, 13 present a qualitative resemblance to a number of equilibrium phase transitions. In Figure 6A, the firing rate $\nu$ vs. $k_r$ of the purple curve ($V_c=-54$mV) exhibits a sharp transition; for $k_r\ge k_r^c=21.15\times 10^{-2}$s^-1 we find power law behavior $(\nu -{\nu }_c)\propto {(k_r-k_r^c)}^{\beta }$ with ${\nu }_c\approx 58$mHz and scaling exponent $\beta \approx 0.32$ (Figure 17). For $V_c>-54$mV the transition appears smoothed out. This resembles the behavior of the magnetization $M$ vs. temperature $T$ for a ferromagnet close to the Curie point. In zero external magnetic field $(H=0)$ the magnetization rises abruptly for $T<T_c$, exhibiting power law behavior $M\propto {(T_c-T)}^{\beta }$; experimentally, the scaling exponent for systems in the Ising universality class is $0.31\le \beta \le 0.33$; for the Ising model in 3D it is $\beta \approx 0.325$ [17].

Figure 17

Figure 17. Scaling of the firing rate as a function of the recovery rate, for Region II of the purple curve in Figure 6A. The slope of the linear fit (blue line) gives a scaling exponent $\beta =0.317$. The critical values $k_r^{(c)}$ and ${\nu }_c$ were determined by starting with values very close to the I $\to$ II transition and making small adjustments to make the points fall into a straight line in the log - log plot.

For finite field $(H\ne 0)$ the transition appears smoothed out in the M - T plane. With the correspondence $\nu \leftrightarrow M$, $k_r\leftrightarrow T$, $V_c\leftrightarrow H$ the plots in Figure 6 resemble the magnetization vs. temperature as the external field is turned on. For the magnetic system, a plot $M$ vs. $H$ would display the phenomenon of hysteresis for $T<T_c$. We therefore ask whether the firing rate $\nu$ vs. clamp voltage $V_c$ could show hysteresis too, for certain values of $k_r$.

To investigate the occurrence of hysteresis in our model of the AA, we simulate the 2D voltage independent system in a slightly different way. We start the system in some initial state with ${\chi }_c$, $k_i$, and $k_r$ fixed, and an initial clamp value $V_c^{(i)}$. The clamp voltage is then increased (“adiabatically”) from this initial value to a final value $V_c^{(f)}$ over a time interval $T$, in uniform increments $(V_c^{(f)}>V_c^{(i)})$. The process is then reversed, with the clamp lowered from $V_c^{(f)}$ to $V_c^{(i)}$ over the time $T$. The firing rate is calculated for each time increment $t$ ($t=T/N$ where $N$ is the number of $V_c$ values sampled between the initial and final values) and plotted as a function of $V_c$ for both the forward and reverse process. Using this protocol, we find that for certain parameter choices there is a difference in firing rate between the forward and reverse processes in the vicinity of the I $\to$ II transition, i.e. a hysteresis loop in the $V_c-\nu$ plane. Figure 18 shows the result for $\chi =0.05$, $k_r=13.4\times 10^{-3}$, $k_i=0.25$. The jump in firing rate corresponds to the I $\to$ II transition, and occurs at a slightly different $V_c$ depending on the direction the system approaches from. Note that the hysteresis loop shown here is not due to the subcritical Hopf bifurcation in the system, which corresponds to the transition II $\to$ III. The existence of hysteresis at a subcritical Hopf bifurcation is well known [16]; here it occurs at slightly larger values of $V_c$ (not shown on Figure 18).

Figure 18

Figure 18. Hysteresis in the firing rate of the simulated 2D system (5), (6). Starting with ${\chi }_c=0.05$, $k_r=13.4\times 10^{-3}$, $k_i=0.25$, and $V_c^{(i)}=-1.718$, the clamp value is increased in increments of $2\times 10^{-6}$ until $V_c^{(f)}=-1.716$ (2,000 total clamp values sampled), staying at each $V_c$ value for $t=20$ so that a firing rate can be calculated. The process is then reversed, with the clamp returning to the initial value through the exact same intermediate values. The plot shows that the precise location of the transition between Region I and II depends on the direction of travel.

3 Experimental results

3.1 Measured rates

The purpose of mapping out the behavior of the AA dynamical system in parameter space is to provide guidance for the experiments in order to realize interesting dynamical behavior and understand the phenomenology. We are interested in knowing which of the dynamical phases can be accessed with the present experimental setup. Further, knowledge of the phase diagram in parameter space will guide the choice of alternative channels to improve the AA. To make progress, we need to know roughly where the present experimental system lies in the phase diagram described in Section 2. Though rates for the KvAP have been measured before [7], these measurements are not easily mapped onto our system, as they are more complex (higher dimensional in parameter space) in order to account for the physical gating properties of the channels; whereas our model only has the minimal complexity needed to retain the core dynamics.The equilibrium opening and closing rates of KvAP in the AA setting have been measured in a previous work [14], so we turn our attention to the inactivation and recovery rates $k_i$ and $k_r$. We measure these rates using a modified version of the system which is voltage clamped in the traditional manner (i.e. $R_c$ = 0). The experiments are then carried out in standard electrophysiological fashion [7, 18], using voltage protocols adapted to obtain the “effective” rates of Figure 1.

To measure $k_i(V)$, the system is first held at the resting voltage $V_r=-120$mV, where almost all channels are in the closed state. At $t=0$ the voltage is stepped up to $V_1=+100$mV and held there for a time $t_1=100$ms. At the end of this time interval most channels are open and only few are inactivated, since at $V_1$ the opening rate is faster than $100$ ms and the inactivation rate considerably slower. At $t=t_1$ the voltage is dropped to a lower (typically negative) value $V_2$, held there for a time $t_2$ $(\sim 1s)$ before being stepped back up to $V_1=+100$ mV. Finally the voltage is returned to the resting state $V_r$ in order to start another measurement. The measured quantity is the clamp current (equal to the channel current if we neglect leak currents). The proportion of open channels at time $t_1$ (immediately before the voltage is stepped to $V_2$) is constant. While the system is held at $V_2$, some channels will inactivate with a rate $k_i(V_2)$, thus the second step to $V_1$ will elicit a smaller current than the first. The ratio of these two current peaks as a function of $V_2$ and $t_2$ allows us to extract the rate $k_i(V)$. Figure 19A shows several current traces which illustrate the protocol. In formulas, the initial state is prepared with $p_i(t=0)\approx 0$ and $p_o(t=0)\approx 1$. Since we want the effective rate $O\to I$ we consider $(d/dt)p_o=-k_i(V)p_o$ for the dynamics while the system is held at $V=V_2$. Therefore after the time $t_2$ we have: $p_o(t_2,V)=p_o(0)\exp [-k_i(V)t_2]$. For a given voltage, the current is $I\propto p_o$ so $I_{\mathit{peak}}/I_0=p_o(t_2,V)/p_o(0)=exp[-k_i(V)t_2]$ where $I_{\mathit{peak}}$ is the peak value of the current when the voltage is stepped to $V_1$ the second time, and $I_0$ the initial peak of the current, when the voltage is stepped to $V_1$ the first time (red trace in Figure 19A). The quantity $-(1/t_2)\ln (I_{\mathit{peak}}/I_0)=k_i(V)$ (where $V=V_2$), obtained from traces as in Figure 19A, is plotted vs. $V$ in Figure 19B, together with a fit to the form $k_i(V)=k_0\exp (\beta V)$ (solid line), from which the parameters $k_0$ and $\beta$ are determined.

Figure 19

Figure 19. (A) Representative time traces of the current corresponding to the voltage clamp protocol used to measure the inactivation rate $k_i(V)$. The initial peak (red trace) gives the maximum current $I_0$, and the ratio of subsequent peaks in comparison gives the ratio of inactivated channels after a time $t_2=1$s spent at the voltage indicated in the legend. $(\mathrm{\Delta }t=t_1+t_2)$ (B) The inactivation rate $k_i(V)$ plotted vs. $V$, obtained from time traces as in (A). The solid line is a fit with an exponential function $k_i(V)=k_0\exp (\beta V)$, returning the values $k_0=0.878$s⁻¹, $\beta =8.13$V⁻¹.

To measure the recovery rate $k_r$, the system is prepared in a state where all channels are inactive, which can be achieved by holding it at a high voltage $V_1=100$mV. The voltage is then stepped down to a $V_2$ below the threshold for firing (between $-120$ to $-80$mV), and held there for a time $\mathrm{\Delta }t$; during this time, a fraction of the channels recover (into the closed state), with rate $k_r(V_2)$; then the voltage is stepped back up to $V_1$. A reference measurement is also taken where the system is held at $V_2$ for a long time ($>20$s) before stepping the voltage back to $V_1$. The ratio of measured current over reference current gives the ratio of open channels (i.e. channels which have recovered from inactivation) as a function of $V_2$ and $\mathrm{\Delta }t$; $k_r(V)$ can be extracted by repeating the process for several values of $V_2$ and $\mathrm{\Delta }t$.

From these measurements, the resulting effective rates for the model of Figure 1 are (in s^-1; $V$ in Volts):

$k_i=0.878\times e^{8.13\times V}{k}_r=0.034\times e^{-11.4\times V}$

These rates can be qualitatively compared to the voltage independent phase diagram of Figure 2A by considering the range of possible voltages of our system. Typical experimental conditions are bounded from below by the resting voltage and above by the Nernst potential: −200 mV $<V<$ 42 mV. This corresponds to a window $0.02<k_r<0.33$ and $0.17<k_i<1.24$ for the rates. This rectangle lies entirely in Region IV of the phase diagram of Figure 2A.Though only heuristic, this comparison is consistent with the absence of autonomous oscillations in the present experimental system. Thus, a single AA with the current setup based on the KvAP channel seems limited to single shot APs. In the next section we report some preliminary measurements on a system of two connected AAs, where we explore the feasibility of obtaining autonomous oscillations.

3.2 Connected AAs

In the future it will be interesting to build networks of interconnected AAs. As a first step, we connected two AAs through electronic “synapses”. Our synapse is a current clamp which takes as command voltage the voltage in the “pre-synaptic” axon and, if this voltage is above a set threshold, delivers a proportional current into the “post-synaptic” axon (similar to a much simplified version of the “dynamic clamp” used in some electrophysiology experiments [19]). Thus a synapse is characterized by two parameters: the threshold voltage $V_T$, and the coefficient of proportionality $\alpha$ between input voltage and output current. A synapse connecting AA1 and AA2 delivers a current into AA2 given by $I_2(t)=\alpha V_1(t)\mathrm{\Theta }[V_1(t)-V_T]$ where $\mathrm{\Theta }$ is the step function. For $\alpha >0$ the synapse is “excitatory” while for $\alpha <0$ it is inhibitory. In the following experiments we keep the threshold at $V_T=0$, while typical values of the synapse “strength” are $|\alpha |\sim 10$pA/mV $=10$nS.

The simplest system that can be made with this setup consists of two AAs connected by a single “excitatory” synapse $(\alpha >0)$. Figure 20 shows experimental time traces of the voltage in the two axons for this configuration. The first axon (AA1: blue line, $V_1$) is caused to fire in the standard way, by stepping up its clamp (CLVC1) from the resting value to a value above threshold; as $V_1$ crosses zero, the synapse starts to inject positive current into AA2, eventually causing it to fire (AA2: red line, $V_2$). As channels inactivate, $V_1$ crosses zero again in the downwards direction, the synapse stops injecting current, and $V_2$ is pulled back to the resting potential by a combination of inactivation and its clamp (CLVC2). During this whole process CLVC2 is held steady at the resting value. $V_1$ does not come back to the resting potential because there are no further inputs to CLVC1 after the initial step up. The end result is that AA2 goes through a complete action potential cycle, including repolarization, with AA1 acting as the input. If AA2 was connected in the same way to a third axon AA3, a similar action potential cycle would be generated in AA3, and so on. A system of several AAs linked in such a way would allow for discrete spatial propagation of action potentials. This configuration is similar to a previously reported result [3] in which the firing of an AP in AA1 caused firing in AA2, the only difference being that AA1 is made to fire via adjustment of its CLVC rather than using an external current source. To summarize: in this configuration, AA1 provides an input signal to AA2, which then fires a complete action potential cycle. We could think of AA1 as a sensory input (which could be realized in practice by embedding light or chemically gated channels in AA1, for example).

Figure 20

Figure 20. Time traces from two AAs connected by one excitatory synapse. Blue is $V_1$ (the voltage in AA1), red is $V_2$. AA1 is caused to fire by raising its CLVC above threshold (to −20 mV); the CLVC of AA2 is held constant at −100 mV. As $V_1$ crosses zero, the synapse starts to inject current into AA2, eventually causing it to fire. As $V_1$ re-crosses zero downwards, the synapse shuts off and $V_2$ is pulled back to the initial resting potential by its CLVC and inactivation. The end result is a complete action potential cycle for AA2, which returns to its initial “resting” state.

As shown in the previous section, as a result of the inactivation dynamics of the KvAP channel, a single AA does not sustain autonomous oscillations for a constant input current. Nevertheless, with two such AAs it is in principle possible to construct an oscillator. For this purpose we add a second synapse to the previous construction, providing inhibitory feedback. Now AA1 connects to AA2 through an excitatory synapse $(\alpha >0)$ and AA2 connects back to AA1 through an inhibitory synapse $(\alpha <0)$. Figure 21A shows corresponding experimental time traces (blue trace: $V_1$; red trace: $V_2$).

Figure 21

Figure 21. (A) Time traces of the voltage from two AAs connected by one excitatory and one inhibitory synapse, in a feedback loop. To initiate the process, CLVC1 is raised above threshold at $t\approx 1.5$s and then held constant; CLVC2 is held fixed at $V_r\approx -100$mV throughout. Firing of AA1 (blue trace) causes a positive synaptic current (${\alpha }_{12}=6$nS, $V_T=0$) which induces firing of AA2 (red trace). When AA2 fires, the negative synaptic current (${\alpha }_{21}=-25$nS, $V_T=0$) injected in AA1 pulls $V_1$ down sharply. In this experiment, not enough channels in AA1 recovered from inactivation (during the negative voltage swing) in order for AA1 to fire again, which would lead to sustained oscillations. (B) Numerical simulation of the system in (A), using the (voltage dependent) model (1), (2), (3) with a combination of measured rates and fitted parameters: $N_0=250$, $C=275$pF, $\chi =167$pS, ${\alpha }_{12}=7.33$nS, ${\alpha }_{21}=-10.67$nS, and $V_T=0$. Note that the discrepancy in synapse strengths is due to a large leak current ${\chi }_{\ell }$ in the experiment.

The rising part is similar to Figure 20, however now when $V_2$ crosses zero, the inhibitory synapse starts injecting negative current into AA1, pulling $V_1$ down to negative values below the resting potential (“hyperpolarization”). As $V_2$ crosses zero again on the falling edge, the synapse shuts off and AA1 repolarizes ($V_1$ rises again) since the clamp CLVC1 is kept steady at the stepped up value (i.e. constant stimulus conditions for AA1). Now in principle the process could repeat and generate a train of APs, i.e. an oscillator. Specifically, if the hyperpolarizing step lasts a sufficient amount of time, depending on synapse strength and channel inactivation/recovery rates, then when the action potential in axon two ceases, the voltage in axon 1 will return to close to the CLVC1 value (which is above threshold) and fire again. This does not quite happen in the experiment shown (but notice the little blip in the traces at $t\approx 3.3$s, which is a partial firing of the system), for the reason that too many channels in AA1 are still inactivated to fire a full AP.

In Figure 21B we show time traces for $V_1$ and $V_2$ from a numerical simulation of the system, using the voltage dependent model (1), (2), (3) for the individual axons (we discuss it further in the next section). We used the experimentally measured rates (see previous section) and fitted the remaining parameters $(N_0,C,{\alpha }_1,{\alpha }_2)$. It is apparent that the model reproduces the dynamics of the real system quite well. We may then interrogate the model on the conditions for obtaining oscillations. It turns out that adjusting synapse strength is enough. Figure 22 shows autonomous oscillations in the model, obtained with the same parameter settings as for Figure 21B, except the strength of the inhibitory synapse AA2 $\to$ AA1 has been increased from ${\alpha }_{21}=-10.67$nS to ${\alpha }_{21}=-20$nS. We conclude that autonomous oscillations are achievable with the current experimental system (though we have not yet been able to obtain them). The key factor in determining whether oscillations are sustained or die out (as in the experiment of Figure 21A) is how low AA1 is pulled by the negative feed back from AA2. A hyperpolarization value for AA1 of $\lesssim -200$ mV is indicated by the simulations to be the minimum required for sustained oscillations under the present conditions.

Figure 22

Figure 22. Time trace from a simulation of the two axon system, showing autonomous oscillations (AP trains). Blue trace is $V_1$, red trace $V_2$. Parameters are the same as for Figure 21B, except the inhibitory synapse strength has been increased to ${\alpha }_{21}=-20$nS. The resulting increased hyperpolarization of AA1 allows more channels to recover from inactivation, resulting in sustained oscillations.

4 Discussion

The dynamical system (1), (2), with voltage dependent rates, describes the dynamics of the physical AA well [4]. The voltage dependence of the rates is taken of the Arrhenius form, $k=\kappa \exp [\alpha (V-V_0)]$. The same model can be used to describe each node of a small network of AAs. This is shown in Section 3.2 where we compare model and experiments for an elaborate experimental system consisting of two AAs connected by “synapses”. Figure 21B shows the time traces obtained numerically for two coupled AAs represented by Equations 1, 2, to be compared to the experimental traces in Figure 21A. The parameter values used were obtained through a combination of measurements and fitting the experimental traces [3, 4]. In Section 2, we discussed dynamics and bifurcations using an even simpler model, with voltage independent inactivation and recovery rates. While the number of parameters is reduced, the simulations show that the phenomenology remains the same. The voltage dependent model produces the same four types of behavior (AP trains, oscillations, damped oscillations, and single shot AP) as we found in the voltage independent model, for a wide range of parameter values. Figure 23 shows example time traces from the system (1), (2), (3) with voltage dependent inactivation and recovery rates: $k_i=3e^{\mathit{20V}}$, $k_r={\kappa }_r{e}^{-20V}$, $V_c=-56$mV, and all other parameters identical to the voltage independent model. In the figure, the system traverses through the four regions in the same fashion as before (I $\to$ II $\to$ III $\to$ IV) as ${\kappa }_r$ is increased; holding ${\kappa }_r$ fixed and varying another parameter again produces the same four regions of behavior. Notably, the sharp transition between Region I and Region II which was found in the voltage independent model is preserved in the full model.

Figure 23

Figure 23. Representative time traces of the model with voltage dependent recovery and inactivation rates in the form of (3). The x-axis is time (s), and the traces from top to bottom are representative of regions I$-$IV, and correspond to ${\kappa }_r=0.037,0.038,0.06$, and $0.3$s^-1, respectively. The fixed inactivation and recovery parameters are: ${\kappa }_i=3$s^-1, ${\alpha }_r=-20$V⁻¹, ${\alpha }_i=20$V⁻¹, $V_0^{(r)}=V_0^{(i)}=0$, with clamp value $V_c=-56$mV. Other parameters ($N_0,C$, etc.) are identical to the voltage independent case (Figure 2A). The top two traces are chosen to show the sharp transition between Region I and Region II.

The simplified Hodgkin-Huxley type model we use captures the dynamics of the experimental system quite well, with a minimum number of parameters. The reduced number of parameters in the model allows us to map out important features of the system’s phase diagram in parameter space (Figure 2A). The main conclusion is that the AA, a synthetic biology system consisting of one voltage gated channel species with inactivation, is dynamically equivalent to the biological system of two voltage gated channel species without inactivation (the Morris-Lecar system). Everything the Morris-Lecar system can do, the AA can do, in principle. This raises the question of whether action potentials dependent on a single gated channel species exist (or have existed) in nature. As far as we know, no such system has been identified thus far.

We have discussed in detail some of the bifurcations which take the system from one region of the phase diagram to another; these have a universal character, which therefore should be maintained across different systems, from the AA to the Morris-Lecar dynamical system to the barnacle muscle fiber to the rat neuron. Indeed, the Hopf bifurcation corresponding to the onset of AP trains, which here we discuss for the AA (Section 2.3), is well established for the neuron [9]. In general, Hopf bifurcations drive the onset of oscillatory behavior in a variety of dynamical systems. One example is protein expression networks [20]; while the mechanism underlying the oscillations varies, it is generally associated with negative feedback loops or time delays [21, 22]. Less established is the transition separating regions I and II in the phase diagram, in fact we do not find that it is discussed in the literature either in theory or experiments. Morris and Lecar mention “small” damped oscillations in their system [8], which may refer to the behavior across regions II and III.

A qualitative analogy with the magnetization transition prompted us to look for hysteresis, which is indeed present.

Hysteresis is an important characteristic of many dynamical systems. It is normally associated with discontinuous transitions; for example, a sequence of transitions in Taylor - Couette flow, associated with discontinuous jumps in the flow pattern, is hysteretic [23]. On the other hand, driven magnetic systems display continuous behavior and hysteresis, due to long relaxation times [24]. Hysteresis may be essential to confer robustness to biological clocks [15, 25].

We have discussed these bifurcations mainly as a function of the inactivation and recovery rates of the channels, $k_i$ and $k_r$, corresponding to the phase diagram of Figure 2A. The reason is to establish guidelines for the future choice of different channels to improve the experimental system. The main conclusion here is that we want a channel with much faster (or more strongly voltage dependent) inactivation. At the same time, it should be noted that the same bifurcations can instead be explored as a function of parameters which are experimentally controlled in the AA. For example, Figure 24 shows plots of the firing rate vs. clamp conductance ${\chi }_c$ for the 2D model of Section 2.3, obtained from the simulation for different values of the clamp voltage $V_c$. The transition Region I $\to$ Region II is again visible as a sharp increase in firing rate as ${\chi }_c$ is lowered past a critical value. This is the same transition as in Figures 6, 13. The fact that this transition is also present in this slice of the phase diagram means that it may be possible to design an experiment to observe this transition in the AA, as both ${\chi }_c$ and $V_c$ are control parameters. It would be difficult to follow this behavior in experiments with living cells, because one does not have the same control over the experimental parameters. Also, a bifurcation which is sharp in the model may appear smoothed out in the more complex environment of the cell. However, the dynamics of the biological system is likely to retain a signature of the underlying bifurcation, which therefore may provide a way to classify the behavior.

Figure 24

Figure 24. Firing rate of the AA as a function of the clamp conductance, ${\chi }_c$, computed from the 2D voltage independent model. Each curves corresponds to a different $V_c$ values (legend). For lower clamp values a sharp transition occurs in firing rate as ${\chi }_c$ is increased past a critical value (Region II $\to$ Region I). Past the transition, the firing rate continues to slowly decrease until a discontinuity occurs due to the clamp current overwhelming the channel current (no firing).

Finally, we present experiments where we connect two AAs through electronic “synapses” - a step towards constructing networks. In the configuration corresponding to Figure 20, a step perturbation of AA1 (the “input”) evokes a single shot AP in AA2 (the “output”). In terms of signal processing, suppose AA1 was a light sensitive channel, and the input consisted of a (slowly) blinking light. The system would encode each event of the light turning on into a single standardized AP (independent of the duration of the “on” phase), which could be further propagated down a network. By connecting two AAs in a feedback loop, it should be possible to construct an autonomous limit cycle oscillator. Figure 21A displays an attempt which was not quite successful, but indicates the feasibility of such a system. The difficulty lies in the fragility of the experimental system, which makes it difficult to tune parameters “on the fly”. Our plan going forward is to both demonstrate this oscillator and also improve robustness.

The parameters which are currently directly controlled in the experiments are the clamp voltage and clamp conductance. Though much more cumbersome, the inactivation and recovery rates for a given channel could in principle be modified as well. Previous work has shown that different compositions of the lipid membrane has an affect on the kinetics of the channel [7], the caveat being that there is no quantitive way of knowing how the rates will change in response to a change in membrane composition. Channel kinetics also change as a function of temperature [26], which could be another way to indirectly tune the system. A combination of such methods may be an effective strategy for exploring the parameter space experimentally. To expand the system beyond a two axon setup, further work is needed due to the difficulty in maintaining multiple functioning AAs, owing to the fragility of the lipid membrane setup. Possible strategies include using hydrogels to stabilize the system [27, 28], or perhaps moving from a suspended lipid system to a supported lipid system.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

ZP: Writing–original draft, Methodology, Investigation, Formal Analysis, Data curation, Conceptualization. GZ: Writing–review and editing, Writing–original draft, Supervision, Formal Analysis, Conceptualization.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. The publication fee is paid by the Department of Physics and Astronomy, UCLA.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2024.1452241/full#supplementary-material

References

1. Gerstner W, Kreiter AK, Markram H, Herz AVM, 94, 12740–1. doi:10.1073/pnas.94.24.127401997).Neural codes: firing rates and beyond, Proc Natl Acad Sci U S A,

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Ariyaratne A, Zocchi G. Toward a minimal artificial axon. J Phys Chem B (2016) 120:6255–63. doi:10.1021/acs.jpcb.6b02578

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Vasquez HG, Zocchi G. Coincidences with the artificial axon. EPL (2017) 119:48003. doi:10.1209/0295-5075/119/48003

CrossRef Full Text | Google Scholar

4. Pi Z, Zocchi G. Critical behavior in the artificial axon. J Phys Commun (2021) 5:125013. doi:10.1088/2399-6528/ac43d0

CrossRef Full Text | Google Scholar

5. Prescott SA, Koninck YD, Sejnowski TJ. Plos Comput Biol (2008) 4:e1000198. doi:10.1371/journal.pcbi.1000198

PubMed Abstract | CrossRef Full Text

6. Hodgkin AL, Huxley AF. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol (Lond.) (1952) 117:500–44. doi:10.1113/jphysiol.1952.sp004764

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Schmidt D, Cross SR, MacKinnon R. A gating model for the archeal voltage-dependent K+ channel KvAP in DPhPC and POPE:POPG decane lipid bilayers. J Mol Biol (2009) 390:902–12. doi:10.1016/j.jmb.2009.05.062

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Morris C, Lecar H. Voltage oscillations in the barnacle giant muscle fiber. Biophys J (1981) 35:193–213. doi:10.1016/s0006-3495(81)84782-0

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Koch C. Biophysics of computation. Oxford University Press (1999).

Google Scholar

10. Ruta V, Jiang Y, Lee A, Chen J, MacKinnon R. Functional analysis of an archaebacterial voltage-dependent K+ channel. Nature (2003) 422:180–5. doi:10.1038/nature01473

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Qi X, Zocchi G. Kink propagation in the artificial axon. EPL (2022) 137:12005. doi:10.1209/0295-5075/ac44e2

CrossRef Full Text | Google Scholar

12. Tsumoto K, Kitajima H, Yoshinaga T, Aihara K, Kawakami H. Bifurcations in morris–lecar neuron model. Neurocomputing (2006) 69:293–316. doi:10.1016/j.neucom.2005.03.006

CrossRef Full Text | Google Scholar

13. Liu C, Liu X, Liu S. Bifurcation analysis of a Morris–Lecar neuron model. Biol Cybernetics (2014) 108:75–84. doi:10.1007/s00422-013-0580-4

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Ariyaratne A, Zocchi G. Nonlinearity of a voltage-gated potassium channel revealed by the mechanical susceptibility. PRX (2013) 3:011010. doi:10.1103/physrevx.3.011010

CrossRef Full Text | Google Scholar

15. Vilar JMG, Kueh HY, Barkai N, Leibler S. Mechanisms of noise-resistance in genetic oscillators. PNAS (2002) 99:5988–92. doi:10.1073/pnas.092133899

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Strogatz SH. Nonlinear dynamics and chaos. Westview Press, Boulder, Colorado: Westview Press (2015).

Google Scholar

17. Goldenfeld N. Lectures on phase transitions and the renormalization group (frontiers in physics v. 85). CRC Press, Boca Raton, Florida: CRC Press (2018).

Google Scholar

18. Hille B. Ion channels of excitable membranes. 3rd ed. (2001) Sinauer, Sunderland, Mass.

Google Scholar

19. Wilders R. Dynamic clamp: a powerful tool in cardiac electrophysiology. The J Physiol (2006) 576:349–59. doi:10.1113/jphysiol.2006.115840

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Xiong LI, Garfinkel A Progress in biophysics and molecular biology, 174 (2022). p. 28.

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Lahav G, Rosenfeld N, Sigal A, Geva-Zatorsky N, Arnold J Levine AJ, Elowitz MB, et al. Dynamics of the p53-Mdm2 feedback loop in individual cells. Nat Genet (2004) 36:147–50. doi:10.1038/ng1293

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Bar-Or RL, Maya R, Segel LA, Alon U, Levine AJ, Oren M. PNAS (2000) 97:11250. doi:10.1073/pnas.210171597

PubMed Abstract | CrossRef Full Text

23. Coles D. Transition in circular Couette flow. J Fluid Mech (1965) 21:385–425. doi:10.1017/s0022112065000241

CrossRef Full Text | Google Scholar

24. Chakrabarti BK, Acharyya M. Dynamic transitions and hysteresis. Rev Mod Phys (1999) 71:847–59. doi:10.1103/revmodphys.71.847

CrossRef Full Text | Google Scholar

25. Barkai N, Leibler S. Circadian clocks limited by noise. Nature (1999) 403:267–8. doi:10.1038/35002258

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Ranjan R, Logette E, Marani M, Herzog M, Tâche V, Scantamburlo E, et al. A kinetic map of the homomeric voltage-gated potassium channel (kv) family. Front Cell Neurosci (2019) 13:358. doi:10.3389/fncel.2019.00358

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Jeon T-J, Malmstadt N, Schmidt JJ. Hydrogel-encapsulated lipid membranes. J Am Chem Soc (2006) 128:42–3. doi:10.1021/ja056901v

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Jeon T-J, Malmstadt N, Poulos JL, Schmidt JJ. Black lipid membranes stabilized through substrate conjugation to a hydrogel. Biointerphases (2008) 3:FA96–FA100. doi:10.1116/1.2948314

PubMed Abstract | CrossRef Full Text | Google Scholar