---
title: LJP Theory and Correction
description: What causes Liquid Junction Potential (LJP) and how to correct for it in patch-clamp experiments
date: 2021-01-13
tags: science
---

# Liquid Junction Potential (LJP) Theory and Correction

This page contains notes about the analytical methods [LJPcalc](https://swharden.com/LJPcalc/) uses to calculate LJP from ion tables
 as well as notes for experimenters about what causes LJP and how to compensate for it in electrophysiology experiments.

> 💡 **Use [LJPcalc](https://swharden.com/LJPcalc/) to calculate Liquid Junction Potential (LJP) in your browser**

## Table of Contents

![](TOC)

## LJP Calculation Notes

### LJPcalc Calculation Method

LJPcalc calculates the liquid junction potential according to the stationary Nernst-Planck equation which is typically regarded as superior to the simpler Henderson equation used by most commercial LJP calculators. Both equations produce nearly identical LJPs, but the Henderson equation becomes inaccurate as ion concentrations increase, and also when calculating LJP for solutions containing polyvalent ions.

### LJPcalc Ion Mobility Library

LJPcalc uses an extensive [ion mobility library](../iontable)

The ion mobility table is stored in Markdown format. Not only does Markdown make it easy to display the table nicely in a browser, 
but it also makes the table easy to edit in any text editor. Users desiring to use their own ion mobilities or add new ions to the table
can do so by editing the `IonTable.md` file adjacent to `LJPcalc.exe` as needed.

### Influence of Ion Sequence on Calculated LJP

> 💡 LJPcalc automatically sorts the ion table into an ideal sequence prior to solving for LJP. Attention only needs to be paid to the ion sequence if automatic sorting is disabled.

When calculating LJP for a set of ions it is important to consider the sequence in which they are listed. Additional information can be found in [Marino et al., 2014](https://arxiv.org/abs/1403.3640) which describes the exact computational methods employed by LJPcalc.

* **The last ion's c0 may be overridden** to achieve electroneutrality on the c0 side. This will not occur if the sum of charge on the c0 side is zero.

* **cL for most ions will be slightly adjusted** to achieve electroneutrality on the cL side. The second-to-last ion's cL (which cannot equal its c0) will remain fixed, while the last cL will be adjusted to achieve electroneutrality. During the solving process all cL values (but the second-from-last) will be slightly adjusted. The adjustments are likely negligible experimentally, but this is why cL values in the output table slightly differ from those given for inputs.

### Effect of Temperature on LJP

**The LJP is temperature dependent.** There are two sources of temperature-dependent variation: the Einstein relation and the conductivity table. The former can be easily defined at calculation time, while the latter requires modifying conductances in the ion mobility table. These modifications typically have a small effect on the LJP, so standard temperature (25C) can be assumed for most applications.

**The [Einstein relation](https://en.wikipedia.org/wiki/Einstein_relation_(kinetic_theory))** defines diffusion as **`D = µ * k * T`** where:

* **`D`** is the diffusion coefficient
* **`µ`** (mu) is [ionic mobility](https://en.wikipedia.org/wiki/Electrical_mobility)
* **`k`** is the [Boltzmann constant](https://en.wikipedia.org/wiki/Boltzmann_constant) (1.380649e-23 J / K)
* **`T`** is temperature (K)

**The ion conductivity table is temperature-specific.** Ion conductivity was measured experimentally and varies with temperature. The ion conductivity table here assumes standard temperature (25C), but ion conductivity values can be found for many ions at nonstandard temperatures. LJPcalc users desiring to perform LJP calculations at nonstandard temperatures are encouraged to build their own temperature-specific ion tables.

### Calculating Ionic Mobility from Charge and Conductivity

Ionic mobility is **`µ = Λ / (N * e² * |z|)`** where:

* **`µ`** (mu) is [ionic mobility](https://en.wikipedia.org/wiki/Electrical_mobility) (m² / V / sec)
* **`Λ`** (Lambda) is [molar conductivity](https://en.wikipedia.org/wiki/Molar_conductivity) (S * cm²/ mol)
* **`N`** is the [Avogadro constant](https://en.wikipedia.org/wiki/Avogadro_constant) (6.02214076e23 particles / mol)
* **`e`** is the [elementary charge](https://en.wikipedia.org/wiki/Elementary_charge) (1.602176634e-19 Coulombs)
* **`z`** is the absolute value of the [elementary charge](https://en.wikipedia.org/wiki/Elementary_charge) of the ion

## How to Correct for LJP in Electrophysiology Experiments

Patch-clamp electrophysiologists impale cells with glass microelectrodes to measure or clamp their voltage.
Amplifier offset voltage is adjusted to achieve a reading of zero volts when the pipette is in open-tip configuration with the bath,
but this voltage includes offset for a _liquid junction potential_ (LJP) caused by the free exchange of ions with different mobilities between the pipette and bath solutions.
Whole-cell patch-clamp experiments typically fill the pipette with large anions like gluconate, aspartate, or methanesulfonate, and their low mobility 
(relative to ions like K, Na, and Cl) causes them to accumulate in the pipette and produce a LJP (where the bath is more positive than then pipette).
After establishment of whole-cell configuration, ions no longer freely move between pipette and bath solutions (they are separated by the cell membrane), so
there is effectively no LJP but the offset voltage is still offsetting as if LJP were present. By knowing the LJP, the scientist can
adjust offset voltage to compensate for it, resulting in more accurate measured and clamped voltages.

<p align="center">V<sub>meter</sub> = V<sub>cell</sub> + LJP</p>

**To correct for LJP,** the electrophysiologist must calculate LJP mathematically (using software like LJPcalc) or estimate it experimentally (see the section on this topic below). Once the LJP is known it can be compensated for experimentally to improve accuracy of recorded and clamped voltages.

<p align="center">V<sub>cell</sub> = V<sub>meter</sub> - LJP</p>

> ⚠️ This method assumes that the amplifier voltage was zeroed at the start of the experiment when the pipette was in open-tip configuration 
with the bath, and that concentration of chloride (if using Ag/AgCl electrodes) in the internal and bath solutions are stable throughout experiments.

<div align="center">

![](ljp.png)

</div>

### Example Patch-Clamp LJP Calculation & Correction

This ion set came from in [Figl et al., 2003](https://medicalsciences.med.unsw.edu.au/sites/default/files/soms/page/ElectroPhysSW/AxoBits39New.pdf) Page 8. They have been loaded into LJPcalc such that the pipette solution is c0 and the bath solution is cL. Note that the order of ions has been adjusted to place the most abundant two ions at the bottom. This is ideal for LJPcalc's analytical method.


 Name       | Charge | pipette (mM) | bath (mM)      
------------|--------|--------------|---------
 K          | +1     | 145          | 2.8
 Na         | +1     | 13           | 145
 Mg         | +2     | 1            | 2
 Ca         | +2     | 0            | 1
 HEPES      | -1     | 5            | 5
 Gluconate  | -1     | 145          | 0           
 Cl         | -1     | 10           | 148.8


Loading this table into LJPcalc produces the following output:

```
Values for cL were adjusted to achieve electro-neutrality:

 Name               | Charge | Conductivity (E-4) | C0 (mM)      | CL (mM)      
--------------------|--------|--------------------|--------------|--------------
 K                  | +1     | 73.5               | 145          | 2.8098265   
 Na                 | +1     | 50.11              | 13           | 144.9794365 
 Mg                 | +2     | 53.06              | 1            | 1.9998212   
 Ca                 | +2     | 59.5               | 0            | 0.9999109   
 HEPES              | -1     | 22.05              | 5            | 4.9990023   
 Gluconate          | -1     | 24.255             | 145          | 0           
 Cl                 | -1     | 76.31              | 10           | 148.789725

Equations were solved in 88.91 ms
LJP at 20 C (293.15 K) = 16.052319631180264 mV
```

> 💡 _[Figl et al., 2003](https://medicalsciences.med.unsw.edu.au/sites/default/files/soms/page/ElectroPhysSW/AxoBits39New.pdf) Page 8 calculated a LJP of 15.6 mV for this ion set (720 µV lesser magnitude than our calculated LJP). As discussed above, differences in ion mobility table values and use of the Nernst-Planck vs. Henderson equation can cause commercial software to report values slightly different than LJPcalc. Experimentally these small differences are negligible, but values produced by LJPcalc are assumed to be more accurate. See [Marino et al., 2014](https://arxiv.org/abs/1403.3640) for discussion._

If we have patch-clamp data that indicates a neuron rests at -48.13 mV, what is its true resting potential? Now that we know the LJP, we can subtract it from our measurement:

<p align="center">V<sub>cell</sub> = V<sub>measured</sub> - LJP</p>

<p align="center">V<sub>cell</sub> = -48.13 - 16.05 mV</p>

<p align="center">V<sub>cell</sub> = -64.18 mV</p>

We now know our cell rests at -64.18 mV.

### Zeroed Voltage = LJP + Two Electrode Half-Cell Potentials

The patch-clamp amplifier is typically zeroed at the start of every experiment when the patch pipette is in open-tip configuration 
with the bath solution. An offset voltage (V<sub>offset</sub>) is applied such that the V<sub>measured</sub> is zero. 
This process incorporates 3 potentials into the offset voltage:

* **liquid junction potential (LJP)** between the **pipette** solution and the **bath** solution (mostly from small mobile ions)
* **half-cell potential (HCP)** between the **reference electrode** and the **bath** solution (mostly from Cl)
* **half-cell potential (HCP)** between the **recording electrode** and the **pipette** solution (mostly from Cl)

When the amplifier is zeroed before to experiments, all 3 voltages are incorporated into the offset voltage. 
Since the LJP is the only one that changes after forming whole-cell configuration with a patched cell (it is eliminated),
it is the only one that needs to be known and compensated for to achieve a true zero offset (the rest remain constant).

However, if the [Cl] of the internal or bath solutions change during the course of an experiment 
(most likely to occur when an Ag/AgCl pellet is immersed in a flowing bath solution), 
the half-cell potentials become significant and affect V<sub>measured</sub> as they change. 
This is why agar bridge references are preferred over Ag/AgCl pellets.
See [Figl et al., 2003](https://medicalsciences.med.unsw.edu.au/sites/default/files/soms/page/ElectroPhysSW/AxoBits39New.pdf) 
for more information about LJPs as they relate to electrophysiological experiments.

### Measuring LJP Experimentally

It is possible to measure LJP experimentally, but this technique is often discouraged because issues with KCl reference electrodes make it difficult 
to accurately measure ([Barry and Diamond, 1970](https://link.springer.com/article/10.1007/BF01868010)). However, experimental measurement 
may be the only option to calculate LJP for solutions containing ions with unknown mobilities.

**To measure LJP Experimentally:**<br>
&bull; **Step 1:** Zero the amplifier with intracellular solution in the bath and in your pipette<br>
&bull; **Step 2:** Replace the bath with extracellular solution<br>
&bull; **Step 3:** The measured voltage is the negative LJP (invert its sign to get LJP)

> **✔️ Confirm no drift is present** by replacing the bath with intracellular solution after step 3 to verify the reading is 0. 
If it is not 0, some type of drift is occurring and the measured LJP is not accurate.

> **❓ Why invert the sign of the LJP?** The LJP measured in step 2 is the LJP of the _pipette_ relative to 
the _bath_, but in electrophysiology experiments convention is to refer to LJP as that of the _bath_ relative to the _pipette_. 
LJPs for experiments using typical ACSF bath and physiological K-gluconate pipette solutions are usually near +15 mV.

> **⚠️ Do not measure LJP using a Ag/AgCl reference electrode!** because mobility will be low when the bath is filled with intracellular solution (physiological intracellular solutions have low [Cl]). 
Use a 3M KCl reference electrode instead, allowing high [K] mobility in intracellular solution and high [Cl] mobility in extracellular solution.

## References

* **[Marino et al. (2014)](https://arxiv.org/abs/1403.3640)** - describes a computational method to calculate LJP according to the stationary Nernst-Planck equation. The JAVA software described in this manuscript is open-source and now on GitHub ([JLJP](https://github.com/swharden/jljp)). Figure 1 directly compares LJP calculated by the Nernst-Planck vs. Henderson equation.

* **[Perram and Stiles (2006)](https://pubs.rsc.org/en/content/articlelanding/2006/cp/b601668e)** - A review of several methods used to calculate liquid junction potential. This manuscript provides excellent context for the history of LJP calculations and describes the advantages and limitations of each.

* **[Shinagawa (1980)](https://www.ncbi.nlm.nih.gov/pubmed/7401663)** _"Invalidity of the Henderson diffusion equation shown by the exact solution of the Nernst-Planck equations"_ - a manuscript which argues that the Henderson equation is inferior to solved Nernst-Planck-Poisson equations due to how it accounts for ion flux in the charged diffusion zone.

* **[Lin (2011)](http://www.sci.osaka-cu.ac.jp/~ohnita/2010/TCLin.pdf)** _"The Poisson The Poisson-Nernst-Planck (PNP) system for ion transport (PNP) system for ion transport"_ - a PowerPoint presentation which reviews mathematical methods to calculate LJP with notes related to its application in measuring voltage across cell membranes.

* **[Nernst-Planck equation](https://en.wikipedia.org/wiki/Nernst%E2%80%93Planck_equation)** (Wikipedia)

* **[Goldman Equation](https://en.wikipedia.org/wiki/Goldman_equation)** (Wikipedia)

* **[EGTA charge and pH](https://www.sciencedirect.com/science/article/pii/S0165027099000369?via%3Dihub#FIG1)** - Empirical determination of EGTA charge state distribution as a function of pH.

* **[LJPCalcWin](https://medicalsciences.med.unsw.edu.au/sites/default/files/soms/page/ElectroPhysSW/JPCalcWin-Demo%20Manual.pdf)** - A Program for Calculating Liquid Junction Potentials

* **[LJP Corrections](http://beenhakkerlab.org/lab-links/ephys/Guides/Theory/Liquid-Junction-Potential-Corrections-Axon.pdf)** (Axon Instruments Application Note) describes how to calculate LJP using ClampEx and LJPCalcWin and also summarizes how to measure LJP experimentally

* **[LJP Corrections](https://medicalsciences.med.unsw.edu.au/sites/default/files/soms/page/ElectroPhysSW/AxoBits39New.pdf)** (Figl et al., AxoBits 39) summarizes LJP and discusses measurement and calculation with ClampEx

---
title: Exploring the Voltage-Clamp Membrane Test
date: 2020-10-11 22:24:00
tags: science, circuit
---

# Exploring the Voltage-Clamp Membrane Test

**By modeling a voltage-clamp amplifier, patch pipette, and cell membrane as a circuit** using free circuit simulation software, I was able to create a virtual patch-clamp electrophysiology workstation and challenge model neurons with advanced voltage-clamp protocols. By modeling neurons with known properties and simulating experimental membrane test protocols, I can write membrane test analysis software and confirm its accuracy by comparing my calculated membrane measurements to the values in the original model. A strong advantage of this method (compared to using physical model cells) is that I can easily change values of any individual component to assess how it affects the accuracy of my analytical methods.

<div class="text-center">

![](whole-cell-voltage-clamp-diagram.png)

</div>

**Instead of modeling a neuron, I modeled the whole patch-clamp system:** the amplifier (with feedback and output filtering), pipette (with an imperfect seal, series resistance, and capacitance), and cell (with membrane resistance, capacitance, and a resting potential). After experimenting with this model for a while I realized that advanced topics (like pipette capacitance compensation, series resistance compensation, and amplifier feedback resistance) become much easier to understand when they are represented as components in a circuit with values that can be adjusted to see how the voltage-clamp trace is affected. Many components of the full model can be eliminated to generate ideal traces, and all models, diagrams, and code shown here can be downloaded from my [membrane test repository](https://github.com/swharden/memtest) on GitHub.

## Circuit Components

### Cell

* **`Vm` (Membrane Potential):** Voltage difference across the neuron's membrane. _Neurons typically maintain a membrane potential near -70 mV. In our model we can simulate this by connecting `Rm` to a -70 mV voltage source instead of grounding it as shown in the diagram above._

* **`Rm` (Membrane Resistance):** The resistance across the cell membrane. _Resistance is inversely correlated with membrane conductivity (influenced primarily by the number of open channels in the membrane). Membrane resistance is sometimes termed "input resistance" because in combination with cell capacitance it determines the time constant of the voltage response to input currents._

* **`Cm` (Membrane Capacitance):** The capacitance of a neuron describes how much charge is required to change its voltage. _Larger cells with more membrane surface area have greater capacitance and require more charge (current times time) to swing their voltage._

* **`Tau` (Membrane Time Constant, τ<sub>cell</sub>):** The membrane time constant describes how fast the cell changes voltage in response to currents across its membrane. This is distinctly different than the voltage clamp time constant which describes how fast the cell changes voltage in response to currents delivered through the patch pipette (dependent on Ra, not Rm). _This metric is best thought of with respect to synaptic currents (not currents delivered through the patch pipette). This is a true biological property of the cell, as it exists even when a pipette is not present to measure it. Membrane time constant is membrane capacitance times membrane resistance. If two cells have the same resistance, the larger one (with greater capacitance) will have a slower membrane time constant._

### Pipette

* **`Ra` (Access Resistance):** The resistance caused by the small open tip of the patch pipette. _If a pipette tip gets clogged this resistance will increase, leading to a failed experiment. Access resistance is the primary contributor to series resistance, but a lesser contributor to input resistance._

* **`Rp` (Pipette Resistance):** Resistance between the amplifier and the tip of the pipette. _Resistance of the solution inside the electrode forms a large component of this resistance, but it is such a low resistance is can often be ignored. Its most important consideration is how it combines with Cp to form a low-pass filter inside the pipette (partially overcome by series resistance compensation) to disproportionately degrade fast voltage-clamp transitions._

* **`Rs` (Seal Resistance):** The resistance formed by the seal between the cell surface and the glass pipette. _Ideal experiments will have high seal resistances in the GΩ range._

* **`Rseries` (Series Resistance):** Sum of all non-biological resistances. Access resistance is the largest contributor to series resistance, but pipette resistance and reference electrode resistance also influences it. _Series resistance is bad for two reasons: it acts as a low-pass filter inside the pipette (reducing magnitude of small transients), and it also acts as a voltage divider in series with membrane resistance (resulting in steady-state voltage error). How impactful each of these are to your experiment is easy to calculate or simulate, and a good experiment will have a membrane / series resistance ratio greater than 10._

* **`TauClamp` (Voltage Clamp Time Constant, τ<sub>clamp</sub>):** The voltage clamp time constant describes how fast the cell changes voltage in response to currents delivered through the patch pipette. _This metric is largely determined by access resistance, and it is typically much smaller than the membrane time constant. It describes the relationship between Ra and Cm, and it does not involve Rm. I consider this measurement purely artificial (not biological) because when a pipette is not in a cell this time constant does not exist._

### Amplifier

* **`Vc` (Command Voltage):** This is the voltage the experimenter tries to move the cell toward. _This isn't always exactly what the cell gets though. First, `Cp` and `Rp` form a small low-pass filter delaying measurement of `Vm`. Similarly, `Ra` and `Cm` form a low-pass filter that delays the clamp system from being able to rapidly swing the voltage of the cell. Finally, `Ra` and `Rm` combine to form a voltage divider, leading the amplifier to believe the cell's voltage is slightly closer to `Vc` than it actually is. Many of these issues can be reduced by capacitance compensation and series resistance compensation._

* **`Vo` (Amplifier Output Voltage):** This voltage exiting the amplifier. It is proportional to the current entering the pipette (passing through Rf according to Ohm's law). _Divide this value by `Rf` to determine the current emitted from the amplifier._

* **`Rf` (Feedback Resistance):** Negative feedback for the amplifier. _The greater the resistance the smaller the noise but the smaller the range of the output. Large resistances >1GΩ are used for single channel recordings and lower resistances <1GΩ are used for whole-cell experiments._

* **`Cf` (Feedback Capacitance):** This capacitor forms an RC low-pass filter with `Rf` to prevent ringing or oscillation. _This is tangentially related to capacitance compensation which uses variable capacitance to a computer-controlled voltage to reduce the effects of `Cp`. The main point of this capacitor here is to stabilize our simulation when `Cp` is added._

* **`Io` (Clamp Current):** Current entering the pipette. _This isn't measured directly, but instead calculated from the amplifier's output voltage (measured by an analog-to-digital converter) and calculated as `Vo/Rf` according to Ohm's law._

## Modeling a Patch-Clamp Experiment in LTSpice

**[LTSpice](https://en.wikipedia.org/wiki/LTspice) is a free analog circuit simulator by Analog Devices.** I enjoy using this program, but only because I'm used to it. For anyone trying to use it for the first time, I'm sorry. Watch a YouTube tutorial to learn how to get up and running with it. [Models used in this project are on GitHub](https://github.com/swharden/memtest) if you wish to simulate them yourself.

<div class="text-center img-border">

![](voltage-clamp-circuit.png)

</div>

**This circuit simulates a voltage clamp membrane test** (square pulses, ±5mV, 50% duty, 20 Hz) delivered through a patch pipette (with no pipette capacitance), a 1GΩ seal, 15 MΩ access resistance, in whole-cell configuration with a neuron resting at -70 mV with 500 MΩ membrane resistance and 150 pF capacitance. The Bessel filter is hooked-up through a unity gain op-amp so it can be optionally probed without affecting the primary amplifier. It's configured to serve as a low-pass filter with a cut-off frequency of 2 kHz.

## Simulating a Membrane Test

**The simulated membrane test shows a typical voltage-clamp trace (green)** which is interesting to compare to the command voltage (red) and the actual voltage inside the cell (blue). Note that although the hardware low-pass filter is connected, the green trace is the current passing through the feedback resistor (Rf). A benefit of this simulation is that we can probe anywhere, and being able to see how the cell's actual voltage differs from the target voltage is enlightening.

<div class="text-center img-border">

![](voltage-clamp-simulation.png)

</div>

**If your clamp voltage does not have sharp transitions**, manually define rise and fall times as non-zero values in the voltage pulse configuration options. Not doing this was a huge trap I fell into. If the rise time and fall time is left at `0`, LTSpice will invent a time for you which defaults to 10%! This slow rise and fall of the clamp voltage pulses was greatly distorting the peaks of my membrane test, impairing calculation of I0, and throwing off my results. When using the PULSE voltage source set the rise and fall times to `1p` (1 picosecond) for ideally sharp edges.

**If saving simulation data consider defining the maximum time step.** Leaving this blank is typically fine for inspecting the circuit within LTSpice, but if you intend to save .raw simulation files and analyze them later with Python (especially when using interpolation to simulate a regular sample rate) define the time step to be a very small number before running the simulation.

## Low-Pass Filtering

**Let's compare the output of the amplifier before and after low-pass filtering.** You can see that the Bessel filter takes the edge off the sharp transient and changes the shape of the curve for several milliseconds. This is an important consideration for analytical procedures which seek to measure the time constant of the decay slope, but I'll leave that discussion for another article.

<div class="text-center img-border">

![](voltage-clamp-filter.png)

</div>

## Calculate Clamp Current from Amplifier Output Voltage

**Patch-clamp systems use a digital-to-analog converter which measures voltage coming out of the amplifier to infer the current being delivered into the pipette.** In other words, the magic ability LTSpice gives us to probe current passing through any resistor in the circuit isn't a thing in real life. Instead, we have to use Ohm's law to calculate it as the ratio of voltage and feedback resistance. 

**Let's calculate the current** flowing into the pipette at the start of this trace when the amplifier's output voltage is -192 mV and our command potential is -75 mV:

```
V = I * R
I = V / R
I = (Vout - Vcmd) / Rf
I = ((-192e-3 V) - (-75e-3 V)) / 500e6 Ω
I = -234 pA
```

Notice I use math to get the difference of `Vout` and `Vcmd`, but in practice this is done at the circuit level using a differential amplifier instead of a unity gain op-amp like I modeled here for simplicity.



## Amplifier Feedback Capacitance

**Let's further explore this circuit by adding pipette capacitance.** I set `Cp` to 100 pF (I know this is a large value) and observed strong oscillation at clamp voltage transitions. This trace shows voltage probed at the output of the Bessel filter.

<div class="text-center img-border">

![](voltage-clamp-Cp100-Cf0.png)

</div>

**A small amount of feedback capacitance reduced this oscillation**. The capacitor `Cf` placed across `Rf` serves as an RC low-pass filter to tame the amplifier's feedback. Applying too much capacitance slows the amplifier's response unacceptably. It was impressive to see how little feedback capacitance was required to change the shape of the curve. In practice parasitic capacitance likely makes design of patch-clamp amplifier headstages very challenging. Experimenting with different values of `Cp` and `Cf` is an interesting experience. Here setting `Cp` to 1 pF largely solves the oscillation issue, but its low-pass property reduces the peaks of the capacitive transients.

<div class="text-center img-border">

![](voltage-clamp-Cp100-Cf1.png)

</div>

## Two-Electrode Giant Squid Axon Model

**I created another model to simulate a giant squid axon studied with a two-electrode system.** It's not particularly useful other than as a thought exercise. By clamping between two different voltages you can measure the difference in current passing through the stimulation resistor to estimate the neuron's membrane resistance. [This model is on GitHub](https://github.com/swharden/memtest) too if you want to change some of the parameters and see how it affects the trace.

<div class="text-center img-border">

![](two-electrode-circuit.png)
![](two-electrode-simulation.png)

</div>

**Let's calculate the squid axon's membrane resistance** from the simulation data just by eyeballing the trace.

```
ΔV = (-65 mV) - (-75 mV) = 10 mV <-- Δ command voltage
ΔI = (5 µA) - (-5 µA) = 10 µA <-- Δ amplifier current
```

```
V = I * R
ΔV = ΔI * Rm
Rm = ΔV / ΔI
Rm = 10e-3 V / 10e-6 A
Rm = 1kΩ <-- calculated membrane resistance
```

## Load LTSpice Simulation Data with Python

**LTSpice simulation data is saved in .raw files can be read analyzed with Python** allowing you to leverage modern tools like numpy, scipy, and matplotlib to further explore the ins and outs of your circuit. I'll discuss membrane test calculations in a future post. Today let's focus on simply getting these data from LTSpice into Python. [Simulation data and full Python code is on GitHub](https://github.com/swharden/memtest). Here we'll analyze the .raw file generated by the whole-cell circuit model above.

```python
# read data from the LTSpice .raw file
import ltspice
l = ltspice.Ltspice("voltage-clamp-simple.raw")
l.parse()

# obtain data by its identifier and scale it as desired
times = l.getTime() * 1e3 # ms
Vcell = l.getData('V(n003)') * 1e3  # mV
Vcommand = l.getData('V(vcmd)') * 1e3  # mV
Iclamp = l.getData('I(Rf)') * 1e12  # pA
```

<div class="text-center">

![](voltage-clamp-simple-fig1.png)

</div>

```python
# plot scaled simulation data
import matplotlib.pyplot as plt

ax1 = plt.subplot(211)
plt.grid(ls='--', alpha=.5)
plt.plot(times, Iclamp, 'r-')
plt.ylabel("Current (pA)")

plt.subplot(212, sharex=ax1)
plt.grid(ls='--', alpha=.5)
plt.plot(times, Vcell, label="Cell")
plt.plot(times, Vcommand, label="Clamp")
plt.ylabel("Potential (mV)")
plt.xlabel("Time (milliseconds)")
plt.legend()

plt.margins(0, .1)
plt.tight_layout()
plt.show()
```

**LTSpice simulation data points are not evenly spaced** in time and may require interpolation to produce data similar to an actual recording which samples data at a regular rate. This topic will be covered in more detail in a later post.

## Membrane Test Analysis

**Let's create an ideal circuit, simulate a membrane test, then analyze the data to see if we can derive original values for access resistance (Ra), cell capacitance (Cm), and membrane resistance (Rm).** I'll eliminate little tweaks like seal resistance, pipette capacitance, and hardware filtering, and proceed with a simple case voltage clamp mode.

<div class="text-center">

![](voltage-clamp-memtest-circuit.png)

</div>

> **⚠️ WARNING:** LTSpice voltage sources have a non-negligible conductance by default, so if you use a voltage source at the base of Rm without a defined resistance you'll have erroneous steady state current readings. Prevent this by defining series resistance to a near infinite value instead of leaving it blank.

Now let's run the simulation and save the output...

<div class="text-center">

![](voltage-clamp-memtest-simulation.png)

</div>

I created a diagram to make it easier to refer to components of the membrane test:

<div class="text-center">

![](voltage-clamp-square-membrane-test.png)

</div>

**Think conceptually about what's happening here:** When the command voltage abruptly changes, `Vcell` and `Vcommand` are very different, so the voltage-clamp amplifier delivers a large amount of current right after this transition. The peak current (`Ipeak`) occurs at time zero relative to the transition. The current change between the previous steady-state current (`Iprev`) and the peak current (`Ipeak`) is only limited by `Ra` (since `Cm` only comes in to play after time passes). Let's call this maximum current change `Id`. With more time the current charges `Cm`, raising the `Vcell` toward (`Vcommand`) at a rate described by `TauClamp`. As `Vcell` approaches `Vcommand` the amplifier delivers less current. Altogether, amplifier current can be approximated by an exponential decay function:

<div class="text-center">

I<sub>t</sub> = I<sub>d</sub> * exp(-t / τ<sub>clamp</sub>) + I<sub>ss</sub>

</div>

### Analyze the Capacitive Transient

**The speed at which `Vcell` changes in response to current delivered through the pipette** is a property of resistance (`Ra`) and capacitance (`Cm`). By studying this curve, we can calculate both. Let's start by isolating one curve. We start by isolating individual capacitive transients:

<div class="text-center">

![](voltage-clamp-simple-fig5.png)

</div>

**Fit each curve to a single exponential function.** I'll gloss over how to do this because it is different for every programming language and analysis software. See my [Exponential Fit with Python](https://swharden.com/blog/2020-09-24-python-exponential-fit/) for details. Basically you'll fit a curve which has 3 parameters: `m`, `tau`, and `b`. You may wish to change the sign of tau depending on the orientation of the curve you are fitting. If your signal is low-pass filtered you may want to fit a portion of the curve avoiding the fastest (most distorted) portion near the peak. If you want to follow along, [code for this project is on GitHub](https://github.com/swharden/memtest).

<div class="text-center">

I<sub>t</sub> = m \* exp(-t / tau) + b


![](mt1.png)

</div>

These are the values I obtained by fitting the curve above:

```
m = 667.070
tau = 2.250
b = -129.996
```

Meaning the curve could be modeled by the equation:
```
I = 667.070 * exp(-t / 2.250) -129.996
```

From these values we can calculate the rest:

* `tau` is one of the fitted parameters and has the same time units as the input data. Don't confuse this value with the cell's time constant (which describes how current across `Rm` changes `Vm`), but instead this value is the time constant of the voltage clamp system (where current across `Ra` changes `Vm`). Because `Ra` is much smaller than `Rm`, this will be a much faster time constant.

* State current (`Iss`) is `b` from the curve fit

* The state current before the step will be called `Iprev`

* Change in old vs. new steady state current will be `Idss`

* Peak current (`Ipeak`) occurs at time zero (when the exponential term is 1) so this is simply `m + b`

* `Id` is peak transient current (difference between `Ipeak` and `Iprev`). Some papers call this `I0`, but other papers use that abbreviation to refer to `Ipeak`, so I'll avoid using that term entirely.

We now have:

```
Iss: -129.996 pA
Iprev: -150.015 pA
Idss: 20.019 pA
Ipeak: 537.074 pA
Id: 687.089 pA
dV: 10 mV
TauClamp: 2.250 ms
```

### Calculate Ra

At time zero, access resistance is the thing limiting our ability to deliver current (`Id`) to a known `ΔV` (10 mV). Therefore we can calculate `Ra` using Ohm's law:

```
V = I * R
ΔV = ΔI * R
R = ΔV / ΔI
Ra = dV / Id
Ra = 10e-3 V / 687.089e-12 A
Ra = 14.554 MΩ <-- pretty close to our model 15 MΩ
```

For now let's call this `Ra`, but note that this is technically `Ra` mixed with a small leakage conductance due to `Rm`. Since `Ra` is so much smaller than `Rm` this small conductance doesn't affect our measurement much. Accuracy of this value will be improved when we apply leak current correction described later on this page.

### Calculate Rm

Now that we know `Ra`, we can revisit the idea that the difference between this steady state current (`Iss`) and the last one (`Iprev`) is limited by the sum of `Rm` and `Ra`. let's use this to calculate `Rm` using Ohm's law:

```
V = I * R
I = V / R
ΔI = ΔV / R
R * ΔI = ΔV
(Ra + Rm) * ΔI = ΔV
Ra * ΔI + Rm * ΔI = ΔV
Rm * ΔI = ΔV - Ra * ΔI
Rm = (ΔV - Ra * ΔI) / ΔI
Rm = (dV - Ra * Idss) / Idss
Rm = (10e-3 V - (14.554e6 Ω * 20.019e-12 A)) / 20.019e-12 A
Rm = 485 MΩ <-- pretty close to our model 500 MΩ
```

Accuracy of this value will be improved when we apply leak current correction described later on this page.

### Calculate Cm from Ra, Rm, and Tau

When we raise the cell's voltage (`Vm`) by delivering current through the pipette (`Ra`), some current escapes through `Rm`. From the cell's perspective when we charge it though, `Ra` and `Rm` are in parallel.

```
tau = R * C
C = tau / R
Cm = tau / (1/(1/Ra + 1/Rm))
Cm = 2.250e-3 sec / (1/(1/14.554e6 Ω + 1/485e6 Ω))
Cm = 159 pF <-- pretty close to our model 150 pF
```

Accuracy of this value will be improved when we apply leak current correction described later on this page.

### Calculate Cm from the Area Under the Curve

**Cell capacitance can alternatively be estimated by measuring the area under the capacitive transient.** This method is frequently used historically, and it is simpler and faster than the method described above because it does not require curve fitting. Each method has its pros and cons (e.g., sensitivity to access resistance, hardware filtering, or resilience in the presence of noise or spontaneous synaptic currents). Rather than compare and contrast the two methods, I'll simply describe the theory underlying how to perform this measurement.

<div class="text-center">

![](mt2.png)

</div>

**After an abrupt voltage transition, all current delivered above the steady state current level goes toward charging the cell,** so by integrating this current over time we can calculate how much charge (`Q`) was delivered. I'll describe this measurement as area under the curve (AUC). When summing these data points yourself be sure to remember to subtract steady state current and divide by the sample rate. [Code for this example is on GitHub](https://github.com/swharden/memtest).

**Charge is measured in Coulombs.** Area under the curve is `1515.412 pA*ms`, but recall that a _femtocoulomb_ is 1pA times 1ms, so it's more reasonable to describe the AUC as `1515.412 fC`. This is the charge required to raise cell's capacitance (`Cm`) by `dV`. The relationship is described by:

```
Q = C * ΔV
C = Q / ΔV
Cm = AUC / ΔV
Cm = 1515.412e-15 C / 10e-3 V
Cm = 1515.412e-15 C / 10e-3 V
Cm = 151.541 pF <-- pretty close to our model 150 pF
```

**This value is pretty close to what we expect**, and I think its accuracy in this case is largely due to the fact that we simulated an ideal unfiltered voltage clamp trace with no noise. Its under-estimation is probably due to the fact that a longer period wasn't used for the integration (which may have been useful for this noise-free simulation, but would not be useful in real-world data). Additional simulation experiments with different combinations of noise and hardware filtering would be an interesting way to determine which methods are most affected by which conditions. Either way, this quick and dirty estimation of whole-cell capacitance did the trick in our model cell.

## Correcting for Leak Current

**Why weren't our measurements exact?** `Rm` leaks a small amount of the `Id` current that passes through `Ra` to charge `Cm`. If you calculate the parallel combined resistance of `Ra` and `Rm` you get `14.56 MΩ` which is pretty much exactly what we measured in our first step and simply called `Ra` at the time. Now that we know the value of both resistances we can calculate a correction factor as the ratio of `Ra` to `Rm` and multiply it by both of our resistances. `Cm` can be corrected by dividing it by the square of this ratio.

```
correction = 1 + Ra / Rm
correction = 1 + 14.554 MΩ / 484.96 MΩ
correction = 1.03

Ra = Ra * correction
Rm = Rm * correction
Cm = Cm / (correction^2)
```

<div class="text-center">

Metric | Model | Measured | Corrected | Error
---|---|---|---|--
Ra|15 MΩ|14.55 MΩ|14.99 MΩ|<1%
Rm|500 MΩ|484.96 MΩ|499.51 MΩ|<1%
Cm (fit)|150 pF|159.20 pF|150.06 pF|<1%

</div>

**This correction is simple** and works well when `Ra/Rm` is small. It's worth noting that an alternative to this correction is to solve for `Ra` and `Rm` simultaneously. The [Membrane Test Algorithms](https://mdc.custhelp.com/app/answers/detail/a_id/17006/~/membrane-test-algorithms) used by pCLAMP calculate `Ra` this way, solving the following equation iteratively using the Newton-Raphson method:

```
Ra^2 - Ra * Rt + Rt * (Tau/Cm) = 0
```

**Overall the values I calculated are within a few percent of expectations,** and I'm satisfied with the calculation strategy summarized here. I am also impressed with what we were able to achieve by modeling a voltage-clamped neuron using a free circuit simulator!

## Use a Voltage-Clamp Ramp to Measure Cm

**It's possible to simulate a voltage-clamp _ramp_ and analyze that trace to accurately measure cell capacitance.** A strong advantage of this method is that it does not depend on `Ra`. Let's start by simulating a 10 mV ramp over 100 ms (50 ms down, 50 ms up). When we simulate this with LTSpice and plot it with Python ([screenshots, data, and code is on GitHub](https://github.com/swharden/memtest)) we find that cell voltage lags slightly behind the clamp voltage.

<div class="text-center">

![](mtramp1.png)

</div>

**During voltage-clamp ramps `Vm` lags behind the command voltage because charging `Cm` is limited by `Ra`.** If we measure the difference in this lag between descending and ascending ramps, we can estimate `Cm` in a way that is insensitive to `Ra`. Stated another way, `Ra` only affects abrupt changes in charging rate. Once the cell is charging at a steady rate, that rate of charge is largely unaffected by `Ra` because the stable charging current is already increased to counteract the previous effect `Ra`. Stated visually, `Ra` only affects the rate of charging at the corners of the V. Therefore, let's proceed ignoring the corners of the V and focus on the middle of each slope where the charging rate is stable (and effect of `Ra` is negligible).

<div class="text-center">

![](mtramp2.png)

</div>

**Analysis is achieved by comparing the falling current to the rising current.** We start separately isolating the falling and rising traces, then reverse one of them and plot the two on top of each other. The left and right edges of this plot represent edges of ramps where the system is still stabilizing to compensate for `Ra`, so let's ignore that part and focus on the middle where the charging rate is stable. We can measure the current lag as half of the mean difference of the two traces. Together with the rate of charge (the rate of the command voltage change) we have everything we need to calculate `Cm`.

<div class="text-center">

![](mtramp3.png)

</div>

```
dI = dQ / dt
dI = Cm * dV / dt
Cm = dI / (dV / dT)
Cm = (59.997e-12 A / 2) / (10e-3 V / 50e-3 sec) <-- 10 mV over 50 ms
Cm = 149.993 pF <-- Our model is 150 pF
```

**This is a fantastic result!** The error we do get is probably the result of a single point of interpolation error while converting the unevenly spaced simulation data to an evenly-spaced array simulating a 20 kHz signal. In this ideal simulation this method of calculating `Cm` appears perfect, but in practice it is highly sensitive to sporadic noise that is not normally distributed (like synaptic currents). If used in the real world each ramp should be repeated many times, and only the quietest sweeps (with the lowest variance in the difference between rising and falling currents) should be used for analysis. However, this is not too inconvenient because this protocol is so fast (10 repetitions per second).

## Summary

**This page described how to model voltage-clamp membrane test sweeps and analyze them to calculate Ra, Cm, and Rm.** We validated our calculations were accurate by matching our calculated values to the ones used to define the simulation. We also explored measuring the area under the curve and using voltage-clamp ramps as alternative methods for determining `Cm`. There are a lot of experiments that could be done to characterize the relationship of noise, hardware filtering, and cell properties on the accuracy of these calculations. For now though, I'm satisfied with what we were able to achieve with free circuit simulation software and basic analysis with Python. [Code for this project is on GitHub](https://github.com/swharden/memtest).

<div class="text-center">

Metric | Model | Calculated | Error
---|---|---|---
Ra|15 MΩ|14.99 MΩ|<1%
Rm|500 MΩ|499.51 MΩ|<1%
Cm (fit)|150 pF|150.06 pF|<1%
Cm (auc)|150 pF|151.541 pF|~1%
Cm (ramp)|150 pF|149.993 pF|<.01%

</div>

## Resources

* [LTSpice models, simulations, and Python code on GitHub](https://github.com/swharden/memtest)

* [Download LTspice](https://www.analog.com/en/design-center/design-tools-and-calculators/ltspice-simulator.html) - a high performance SPICE simulation software, schematic capture, and waveform viewer with enhancements and models for easing the simulation of analog circuits.

* The [LC Filter Design Tool](https://rf-tools.com/lc-filter/) makes it easy to design filter circuits using common component values

* [The Patch-clamp Technique Explained And Exercised With The Use Of Simple Electrical Equivalent Circuits](https://mdc.custhelp.com/euf/assets/images/KB864_ModelCells.pdf) by Dirk L Ypey and Louis J. DeFelice
 
* [Series Resistance Compensation](Drexel_Gao_Lab_Series_Resistance_Compensation.pdf) by [Wen-Jun Gao](https://drexel.edu/medicine/about/departments/neurobiology-anatomy/research/gao-lab/)

* [How to correct for series resistance (and whole cell capacitance) in real cells](http://www.billconnelly.net/?p=616) by [Bill Connelly](https://www.utas.edu.au/profiles/staff/health/bill-connelly)

* [Series Resistance. Why it’s bad.](http://www.billconnelly.net/?p=310) by [Bill Connelly](https://www.utas.edu.au/profiles/staff/health/bill-connelly)

* [MultiClamp 700B Theory and Operation](https://mdc.custhelp.com/euf/assets/content/MultiClamp_700B_Manual2.pdf)

* [Introduction to Operational Amplifiers with LTSpice](https://learn.sparkfun.com/tutorials/introduction-to-operational-amplifiers-with-ltspice/all)

* [pyABF](https://swharden.com/pyabf/) - A simple Python interface for ABF files

* [What we talk about when we talk about capacitance measured with the voltage-clamp step method](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3273682/pdf/10827_2011_Article_346.pdf) by Adam L. Taylor (2012)

* [Membrane Capacitance Measurements Revisited: Dependence of Capacitance Value on Measurement Method in Nonisopotential Neurons](https://journals.physiology.org/doi/pdf/10.1152/jn.00160.2009) (Golowasch et al., 2009) discusses 3 ways to calculate capacitance: a current-clamp step, a voltage-clamp step, and a V-shaped pair or voltage-clamp ramps.

* [Techniques for Membrane Capacitance Measurements](https://link.springer.com/chapter/10.1007/978-1-4419-1229-9_7) (Single-Channel Recording, chapter 7) describes the membrane test using a simplified circuit similar to what is discussed here. This is the closest text to a step-by-step guide I've found analyzing the traditional voltage-clamp step protocol.

* [Letter to the editor: Accurate cell capacitance determination from a single voltage step: a reminder to avoid unnecessary pitfalls](https://journals.physiology.org/doi/pdf/10.1152/ajpheart.00503.2016) (Platzer and Zorn-Pauly) discusses the advantages of properly fitting the voltage-clamp curve and calculating I0 instead of just measuring the peak or taking the area under the curve to be the charge.

* [Membrane Test Guide](https://mdc.custhelp.com/app/answers/detail/a_id/17005/~/membrane-test-guide) by Axon / Molecular Devices

* [Membrane Test Algorithms](https://mdc.custhelp.com/app/answers/detail/a_id/17006/~/membrane-test-algorithms) by Axon / Molecular Devices

* [Electronics for electrophysiologists](https://www.biologie.ens.fr/perso/barbour/electronics_for_electrophysiologists.pdf) by Boris Barbour

---
title: Microcontroller Action Potential Generator
date: 2017-08-20 15:12:33
tags: science, circuit, microcontroller, python
---

# Microcontroller Action Potential Generator

__Here I demonstrate how to use a _single_ microcontroller pin to generate action-potential-like waveforms. __The output is similar [my fully analog action potential generator circuit](https://www.swharden.com/wp/2017-08-12-analog-action-potential-generator-circuit/), but the waveform here is created in an entirely different way. A microcontroller is at the core of this project and determines when to fire action potentials. Taking advantage of the pseudo-random number generator ([rand() in AVR-GCC's stdlib.h](http://www.nongnu.org/avr-libc/user-manual/group__avr__stdlib.html#gae23144bcbb8e3742b00eb687c36654d1)), I am able to easily produce unevenly-spaced action potentials which more accurately reflect those observed in nature. This circuit has a potentiometer to adjust the action potential frequency (probability) and another to adjust the amount of overshoot (afterhyperpolarization, AHP). I created this project because I wanted to practice designing various types of action potential _measurement_ circuits, so creating an action potential _generating_ circuit was an obvious perquisite.

![](https://www.youtube.com/embed/2s8t3UsONFs)

__The core of this circuit is a capacitor which is charged and discharged by toggling a microcontroller pin between high, low, and high-Z states.__ In the high state (pin configured as output, clamped at 5V) the capacitor charges through a series resistor as the pin sources current. In the low state (pin configured as output, clamped at 0V) the capacitor discharges through a series resistor as the pin sinks current. In the high-Z / high impedance state (pin configured as an _input_ and little current flows through it), the capacitor rests. By spending most of the time in high-Z then rapidly cycling through high/low states, triangular waveforms can be created with rapid rise/fall times. Amplifying this transient and applying a low-pass filter using a single operational amplifier stage of an [LM-358](http://www.ti.com/lit/ds/symlink/lm158-n.pdf) shapes this transient into something which resembles an action potential. Wikipedia has a section describing how to [use an op-amp to design an active low-pass filter](https://en.wikipedia.org/wiki/Low-pass_filter#Active_electronic_realization) like the one used here.

<div class="text-center img-border">

[![](action-potential-generator-circuit_thumb.jpg)](action-potential-generator-circuit.jpg)

</div>

__The code to generate the digital waveform is very straightforward.__ I'm using PB4 to charge/discharge the capacitor, so the code which actually fires an action potential is as follows:

```c
// rising part = charging the capacitor
DDRB|=(1<<PB4); // make output (low Z)
PORTB|=(1<<PB4); // make high (5v, source current)
_delay_ms(2); // 2ms rise time

// falling part
DDRB|=(1<<PB4); // make output (low Z)
PORTB&=~(1<<PB4); // make low (0V, sink current)
_delay_ms(2); // 2ms fall time
_delay_us(150); // extra fall time for AHP

// return to rest state
DDRB&=~(1<<PB4); // make input (high Z)
```

__Programming the microcontroller__ was accomplished after it was soldered into the device using test clips attached to my ICSP ([USBtinyISP](https://www.ebay.com/sch/i.html?&_nkw=USBtinyISP)). I only recently started using test clips, and for one-off projects like this it's so much easier than adding header sockets or even wiring up header pins.

<div class="text-center img-border">

[![](ap-generator-programmer_thumb.jpg)](ap-generator-programmer.jpg)
[![](ap-generator-programmer-close_thumb.jpg)](ap-generator-programmer-close.jpg)

</div>

__I am very pleased with how well this project turned out!__ I now have an easy way to make irregularly-spaced action potentials, and have a great starting point for future projects aimed at _measuring_ action potential features using analog circuitry.


<div class="text-center img-border">

[![](ap-generator-running_thumb.jpg)](ap-generator-running.jpg)
[![](ap-generator-running-2_thumb.jpg)](ap-generator-running-2.jpg)

</div>

### Notes

*   Action potential half-width (relating to the speed of the action potential) could be adjusted _in software_ by reducing the time to charge and discharge the capacitor. A user control was not built in to the circuit shown here, however it would be very easy to allow a user to switch between regular and fast-spiking action potential waveforms.
*   I am happy that using the 1n4148 diode on the positive input of the op-amp works, but using two 100k resistors (forming a voltage divider floating around 2.5V) at the input and reducing the gain of this stage may have produced a more reliable result.
*   Action potential frequency (probability) is currently detected by sensing the analog voltage output by a rail-to-rail potentiometer. However, if you sensed a noisy line (simulating random excitatory and inhibitory synaptic input), you could easily make an integrate-and-fire model neuron which fires in response to excitatory input.
*   Discussion related to the nature of this "model neuron" with respect to other models (i.e., Hodgkin–Huxley) are on the [previous post](https://www.swharden.com/wp/2017-08-12-analog-action-potential-generator-circuit/).
*   Something like this would make an interesting science fair project

### Source Code on GitHub

*   <https://github.com/swharden/AVR-projects/>

*   <https://github.com/swharden/AVR-projects/tree/master/ATTiny85%202017-08-19%20action%20potential%20generator>

---
title: Action Potential Generator Circuit
date: 2017-08-12 23:04:32
tags: science, circuit, microcontroller, python
---

# Action Potential Generator Circuit

__Few biological cells are as interesting to the electrical engineer as the neuron. Neurons are essentially capacitors (with a dielectric cell membrane separating conductive fluid on each side) with parallel charge pumps, leak currents, and nonlinear voltage-dependent currents. __When massively parallelized, these individual functional electrical units yield complex behavior and underlie consciousness. The study of the electrical properties of neurons ([neurophysiologically](https://en.wikipedia.org/wiki/Neurophysiology), a subset of [electrophysiology](https://en.wikipedia.org/wiki/Electrophysiology)) often involves the development and use of sensitive electrical equipment aimed at studying these small potentials produced by neurons and currents which travel through channels embedded in their membranes. It seems neurophysiology has gained an emerging interest from the hacker community, as evidenced by the success of [Back Yard Brains](https://backyardbrains.com/), projects like the [OpenEEG](http://openeeg.sourceforge.net/doc/hw/), and Hack-A-Day's recent feature [_The Neuron - a Hacker's Perspective_](http://hackaday.com/2017/06/02/the-neuron-a-hackers-perspective/).

While contemplating designs for action potential detection and analysis circuitry, I realized that it would be beneficial to be able to _generate_ action-potential-like waveforms on my workbench. The circuit I came up with to do this is a fully analog (technically mixed signal) action potential generator which produces lifelike action potentials.

![](https://www.youtube.com/embed/MBwQrIVJqfs)

__Cellular Neurophysiology for Electrical Engineers (in 2 sentences):__ Neuron action potentials (self-propagating voltage-triggered depolarizations) in individual neurons are measured in scientific environments using single cell recording tools such as [sharp microelectrodes](http://www.scholarpedia.org/article/Intracellular_recording) and [patch-clamp pipettes](https://en.wikipedia.org/wiki/Patch_clamp). Neurons typically rest around -70mV and when depolarized (typically by external excitatory input) above a threshold they engage in a self-propagating depolarization until they reach approximately +40mV, at which time a self-propagating repolarization occurs (often over-shooting the initial rest potential by several mV), then the cell slowly returns to the rest voltage so after about 50ms the neuron is prepared to fire another action potential. Impassioned budding electrophysiologists may enjoy further reading _[Active Behavior of the Cell Membrane](http://www.bem.fi/book/04/04.htm) _and [_Introduction to Computational Neuroscience_.](http://ecee.colorado.edu/~ecen4831/cnsweb/cns1.html)

<div class="text-center img-border">

![](a.gif)

</div>

__The circuit I describe here produces waveforms which visually mimic action potentials__ rather than serve to replicate the exact conductances real neurons employ to exhibit their complex behavior. It is worth noting that numerous scientists and engineers have designed more physiological electrical representations of neuronal circuitry using discrete components. In fact, the [_Hodgkin-Huxley model_](https://en.wikipedia.org/wiki/Hodgkin%E2%80%93Huxley_model) of the initiation and propagation of action potentials earned Alan Hodgkin and Andrew Huxley the Nobel Prize in Physiology and Medicine in 1936. Some resources on the internet describe how to design lifelike action potential generating circuits by mimicking the endogenous ionic conductances which underlie them, notably _[Analog and Digital Hardware Neural Models](https://people.ece.cornell.edu/land/PROJECTS/NeuralModels/)_, [_Active Cell Model_](https://courses.cit.cornell.edu/bionb442/labs/f2007/lab6.html), and [_Neuromorphic Silicon Neuron Circuits_](http://journal.frontiersin.org/article/10.3389/fnins.2011.00073/full). My goal for this project is to create waveforms which _resemble_ action potentials, rather than waveforms which truly _model_ them. I suspect it is highly unlikely I will earn a Nobel Prize for the work presented here.

<div class="text-center img-border">

[![](circuit_hand_thumb.jpg)](circuit_hand.jpg)

</div>

__The analog action potential simulator circuit I came up with creates a continuous series action potentials.__ This is achieved using a [555 timer](https://en.wikipedia.org/wiki/555_timer_IC) (specifically the [NE555](http://www.ti.com/lit/ds/slfs022i/slfs022i.pdf)) in an astable configuration to provide continuous square waves (about 6 Hz at about 50% duty). The rising edge of each square wave is isolated with a diode and used to charge a capacitor*. While the charge on the capacitor is above a certain voltage, an [NPN transistor](https://en.wikipedia.org/wiki/Bipolar_junction_transistor) (the [2N3904](https://www.onsemi.com/pub/Collateral/2N3903-D.PDF)) allows current to flow, amplifying this transient input current. The capacitor* discharges predictably ([as an RC circuit](https://en.wikipedia.org/wiki/RC_circuit)) through a leak resistor. A large value leak resistor slows the discharge and allows that signal's transistor to flow current for a longer duration. By having two signals (fast and slow) using RC circuits with different resistances (smaller and larger), the transistors are on for different durations (shorter and longer). By making the short pulse positive (using the NPN in [common collector configuration](https://en.wikipedia.org/wiki/Common_collector)) and the longer pulse negative (using the NPN in [common emitter configuration](https://en.wikipedia.org/wiki/Common_emitter)), a [resistor voltage divider](https://en.wikipedia.org/wiki/Voltage_divider) can be designed to scale and combine these signals into an output waveform a few hundred mV in size with a 5V power supply. Pictured below is the output of this circuit realized on a breadboard. The blue trace is the output of the 555 timer.

<div class="text-center img-border">

[![](scope_b_thumb.jpg)](scope_b.jpg)
[![](scope_a_thumb.jpg)](scope_a.jpg)

</div>

Between the capacitance of the rectification diode, input capacitance of the transistor, and stray parasitic capacitance from the physical construction of my wires and the rails on my breadboard, there is sufficient capacitance to accumulate charge which can be modified by changing the value of the leak resistor.


<div class="text-center img-border">

[![](breadboard_thumb.jpg)](breadboard.jpg)

</div>

__This circuit produces similar output when simulated.__ I'm using [LTspice](http://www.linear.com/designtools/software/) (free) to simulate this circuit. The circuit shown is identical to the one hand-drawn and built on the breadboard, with the exception that an additional 0.1 µF capacitor to ground is used on the output to smooth the signal. On the breadboard this capacitance-based low-pass filtering already exists due to the capacitive nature of the components, wires, and rails.

<div class="text-center img-border">

[![](circuit_thumb.jpg)](circuit.png)
[![](sim_trace2_thumb.jpg)](sim_trace2.png)
[![](sim_trace_thumb.jpg)](sim_trace.png)

</div>

__A few improvements naturally come to mind__ when considering this completed, functional circuit:

* **Action potential frequency:** The resistor/capacitor network on the 555 timer determines the rate of square pulses which trigger action potentials. Changing these values will cause a different rate of action potential firing, but I haven't attempted to push it too fast and suspect the result would not be stable is the capacitors are not given time to fully discharge before re-initiating subsequent action potentials.

* **Microcontroller-triggered action potentials:** Since action potentials are triggered by any 5V rising edge signal, it is trivially easy to create action potentials from microcontrollers! You could create some very complex firing patterns, or even "reactive" firing patterns which respond to inputs. 

For example, add a <a href="https://www.adafruit.com/product/439">TSL2561 I2C digital light sensor</a> and you can have a light-to-frequency action potential generator!

* **Adjusting size and shape of action potentials:** Since the waveform is the combination of two waveforms, you can really only adjust the duration (width) or amplitude (height) of each individual waveform, as well as the relative proportion of each used in creating the summation. Widths are adjusted by changing the leak resistor on the base of each transistor, or by adding additional capacitance. Amplitude and the ratio of each signal may be adjusted by changing the ratio of resistors on the output resistor divider.

* **Producing -70 mV (physiological) output:** The current output is electirically decoupled (through a series capacitor) so it can float at whatever voltage you bias it to. Therefore, it is easy to "pull" in either direction. Adding a 10k potentiometer to bias the output is an easy way to let you set the voltage. A second potentiometer gating the magnitude of the output signal will let you adjust the height of the output waveform as desired.

* **The 555 could be replaced by an inverted ramp (sawtooth):** An inverted ramp / sawtooth pattern which produces rapid 5V rising edges would drive this circuit equally well. A <a href="http://www.learningaboutelectronics.com/Articles/Ramp-generator-circuit-with-transistors.php">fully analog ramp generator circuit</a> can be realized with 3 transistors: essentially a constant current capacitor charger with a threshold-detecting PNP/NPN discharge component.

* **This action potential is not all-or-nothing:** In real life, small excitatory inputs which fail to reach the action potential threshold do not produce an action potential voltage waveform. This circuit uses 5V rising edges to produce action potential waveforms. However, feeding a 1V rising edge would produce an action potential 1/5 the size. This is not a physiological effect. However, it is unlikely (if not impossible) for many digital signal sources (i.e., common microcontrollers) to output anything other than sharp rising edge square waves of fixed voltages, so this is not a concern for my application.

* **Random action potentials:** When pondering how to create randomly timed action potentials, the issue of how to generate random numbers arises. This is surprisingly difficult, especially in embedded devices. If a microcontroller is already being used, consider <a href="http://makezine.com/projects/really-really-random-number-generator/">Make's write-up on the subject</a>, and I think personally I would go with a <a href="http://holdenc.altervista.org/avalanche/">transistor-based avalanche nosie generator</a> to create the randomness.

* **A major limitation is that irregularly spaced action potentials have slightly different amplitudes.**I found this out the next day when I created a hardware random number generator (yes, that happened) to cause it to fire regularly, missing approximately half of the action potentials. When this happens, breaks in time result in a larger subsequent action potential. There are several ways to get around this, but it's worth noting that the circuit shown here is best operated around 6 Hz with only continuous regularly-spaced action potentials.

<div class="text-center img-border img-medium">

[![](File_000-2_thumb.jpg)](File_000-2.jpeg)

</div>

__In the video I also demonstrate how to record the output of this circuit using a high-speed (44.1 kHz) 16-bit analog-to-digital converter you already have (the microphone input of your sound card).__ I won't go into all the details here, but below is the code to read data from a WAV file and plot it as if it were a real neuron. The graph below is an _actual recording_ of the circuit described here using the microphone jack of my sound card.

```python
import numpy as np
import matplotlib.pyplot as plt
Ys = np.memmap("recording.wav", dtype='h', mode='r')[1000:40000]
Ys = np.array(Ys)/max(Ys)*150-70
Xs = np.arange(len(Ys))/44100*1000
plt.figure(figsize=(6,3))
plt.grid(alpha=.5,ls=':')
plt.plot(Xs,Ys)
plt.margins(0,.1)
plt.title("Action Potential Circuit Output")
plt.ylabel("potential (mV)")
plt.xlabel("time (ms)")
plt.tight_layout()
plt.savefig("graph.png")
#plt.show()

```

<div class="text-center img-medium">

[![](graph_thumb.jpg)](graph.png)

</div>

__Let's make some noise!__ Just to see what it would look like, I created a circuit to generate slowly drifting random noise. I found this was a non-trivial task to achieve in hardware. Most noise generation circuits create random signals on the RF scale (white noise) which when low-pass filtered rapidly approach zero. I wanted something which would slowly drift up and down on a time scale of seconds. I achieved this by creating 4-bit pseudo-random numbers with a shift register ([74HC595](https://www.sparkfun.com/datasheets/IC/SN74HC595.pdf)) clocked at a relatively slow speed (about 200 Hz) having essentially random values on its input. I used a [74HC14](https://assets.nexperia.com/documents/data-sheet/74HC_HCT14.pdf) inverting buffer (with Schmidt trigger inputs) to create the low frequency clock signal (about 200 Hz) and an extremely fast and intentionally unstable square wave (about 30 MHz) which was sampled by the shift register to generate the "random" data. The schematic illustrates these points, but note that I accidentally labeled the 74HC14 as a 74HC240. While also an inverting buffer the [74HC240](https://assets.nexperia.com/documents/data-sheet/74HC_HCT240.pdf) will not serve as a good RC oscillator buffer because it does not have Schmidt trigger inputs.

An inverting buffer created a fast and a slow clock to produce 4-bit pseudo-random numbers:
<div class="text-center img-border img-medium">

[![](File_000-5_thumb.jpg)](File_000-5.jpeg)

</div>

reminder how an inverting buffer can act as an oscillator:
<div class="text-center">

[![](schmitt-trigger-osc_thumb.jpg)](schmitt-trigger-osc.png)

</div>

the full circuit realized on the breadboard:
<div class="text-center img-border">

[![](File_004_thumb.jpg)](File_004.jpeg)

</div>

output of the 4-bit pseudo-random number generator:
<div class="text-center img-border">

[![](File_001_thumb.jpg)](File_001.jpeg)

</div>

4-bit output smoothed through a single-stage RC filter:
<div class="text-center img-border">

[![](File_002_thumb.jpg)](File_002.jpeg)

</div>

noise combined with action potential waveforms:
<div class="text-center img-border">

[![](File_003_thumb.jpg)](File_003.jpeg)

</div>

__The addition of noise was a success, from an electrical and technical sense.__ It isn't particularly physiological. Neurons would fire differently based on their resting membrane potential, and the peaks of action potential should all be about the same height regardless of the resting potential. _However_ if one were performing an electrical recording through a patch-clamp pipette in perforated patch configuration (with high resistance between the electrode and the internal of the cell), a sharp microelectrode (with high resistance due to the small size of the tip opening), or were using electrical equipment or physical equipment with amplifier limitations, one could imagine that capacitance in the recording system would overcome the rapid swings in cellular potential and result in "noisy" recordings similar to those pictured above. They're not physiological, but perhaps they're a good electrical model of what it's like trying to measure a physiological voltage in a messy and difficult to control experimental environment.

__This project was an interesting exercise in analog land,__ and is completed sufficiently to allow me to move toward my initial goal: creating advanced action potential detection and measurement circuitry. There are many tweaks which may improve this circuit, but as it is good enough for my needs I am happy to leave it right where it is. If you decide to build a similar circuit (or a vastly different circuit to serve a similar purpose), send me an email! I'd love to see what you came up with.

## UPDATE: add a microcontroller

I enhanced this project by creating a [microcontroller controlled action potential generator](https://www.swharden.com/wp/2017-08-20-microcontroller-action-potential-generator/). That article is here: <https://www.swharden.com/wp/2017-08-20-microcontroller-action-potential-generator/>