The personal website of Scott W Harden
Source code for Exploring the Membrane Test with a Voltage-Clamped Neuron Model
---
title: Exploring the Membrane Test with a Voltage-Clamped Neuron Model
date: 2020-10-11 22:24:00
tags:
  - science
  - circuit
---

# Exploring the Membrane Test with a Voltage-Clamped Neuron Model

**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="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="center 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="center 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="center 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="center 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="center 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="center 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="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="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="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="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="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="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="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="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="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="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="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="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="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 ](file:///C:/Users/scott/Downloads/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
Source code for The New Age of QRSS
---
title: The New Age of QRSS
date: 2020-10-03 22:08:00
tags: 
  - qrss
  - amateur radio
---

# The New Age of QRSS

**QRSS is an experimental radio mode that uses frequency-shift-keyed (FSK) continuous wave (CW) Morse code to transmit messages that can be decoded visually by inspecting the radio frequency spectrogram.** The name "QRSS" is a derivation of the [Q code](https://en.wikipedia.org/wiki/Q_code) "QRS", a phrase Morse code operators send to indicate the transmitter needs to slow down. The extra "S" means slow way, _way_ down, and at the typical speed of 6 second dots and 18 second dashes most QRSS operators have just enough time to send their call sign once every ten minutes (as required by federal law). These slow Morse code messages can be decoded by visual inspection of spectrograms created by computer software processing the received audio. A QRSS grabber is a radio/computer setup configured to upload the latest radio spectrogram to the internet every 10 minutes. [QRSS Plus](https://swharden.com/qrss/plus/) is an automatically-updating list of active QRSS grabbers around the world, allowing the QRSS community to see QRSS transmitters being detected all over the world. 

> **TLDR: Get Started with QRSS**
> * Tune your radio to 10.140 MHz (10.1387 MHz USB)
> * Install spectrogram software like [FSKview](https://swharden.com/software/FSKview)
> * Inspect the spectrogram to decode callsigns visually
> * Join the [QRSS Knights](https://groups.io/g/qrssknights) mailing list to learn what's new
> * Go to [QRSS Plus](https://swharden.com/qrss/plus/) to see QRSS signals around the world
> * Design and build a circuit (or [buy a kit](https://www.qrp-labs.com/)) to transmit QRSS

## What is QRSS?

**QRSS allows miniscule amounts of power to send messages enormous distances.** For example, 200 mW QRSS transmitters are routinely spotted on QRSS grabbers thousands of miles away. The key to this resilience lies in the fact that spectrograms can be designed which average several seconds of audio into each pixel. By averaging audio in this way, the level of the noise (which is random and averages toward zero) falls below the level of the signal, allowing visualization of signals on the spectrogram which are too deep in the noise to be heard by ear.

<div class="center border">

![](lopora-signals2.jpg)

</div>

**If you have a radio and a computer, you can view QRSS!** Connect your radio to your computer's microphone, then run a spectrogram like [FSKview](https://swharden.com/software/FSKview/) to visualize that audio as a spectrogram. The most QRSS activity is on 30m within 100 Hz of 10.140 MHz, so set your radio to upper sideband (USB) mode and tune to 10.1387 MHz so QRSS audio will be captured as 1.4 kHz audio tones. 

**[FSKview](https://swharden.com/software/FSKview/) is radio frequency spectrogram software for viewing QRSS and WSPR simultaneously.** I wrote [FSKview](https://swharden.com/software/FSKview/) to be simple and easy to use, but it's worth noting that [Spectrum Lab](https://www.qsl.net/dl4yhf/spectra1.html),  [Argo](https://digilander.libero.it/i2phd/argo/), [LOPORA](https://www.qsl.net/pa2ohh/11lop.htm), and [QRSSpig](https://gitlab.com/hb9fxx/qrsspig) are also popular spectrogram software projects used for QRSS, with the last two supporting Linux and suitable for use on the Raspberry Pi.

<div class="center">

![](fskview.png)

</div>

## QRSS Transmitter Design

**QRSS transmitters can be extraordinarily simple because they just transmit a single tone which shifts between two frequencies.** The simplicity of QRSS transmitters makes them easy to assemble as a kits, or inexpensively designed and built by those first learning about RF circuit design. The simplest designs use a crystal oscillator (typically a [Colpitts](https://en.wikipedia.org/wiki/Colpitts_oscillator) configuration) followed by a buffer stage and a final amplifier (often [Class C](https://en.wikipedia.org/wiki/Power_amplifier_classes#Class_C) configuration using a 2N7000 N-channel MOSFET or 2N2222 NPN transistor). Manual frequency adjustments are achieved using a variable capacitor, supplemented in this case with twisted wire to act as a simple but effective variable capacitor for fine frequency tuning within the 100 Hz QRSS band. Frequency shift keying to transmit call signs is typically achieved using a microcontroller to adjust voltage on a reverse-biased diode (acting as a [varactor](https://en.wikipedia.org/wiki/Varicap)) to modulate capacitance and shift resonant frequency of the oscillator. Following a low-pass filter (typically a 3-pole [Chebyshev](https://en.wikipedia.org/wiki/Chebyshev_filter) design) the signal is then sent to an antenna.

<div class="center medium">

![](qrp-labs-kit-schematic.jpg)

</div>

**[QRP Labs](https://www.qrp-labs.com/) is a great source for QRSS kits.** The kit pictured above and below is one of their earliest kits (the 30/40/80/160m QRSS Kit), but they have created many impressive new products in the last several years. Some of their more advanced QRSS kits leverage things like direct digital synthesis (DDS), GPS time synchronization, and the ability to transmit additional digital modes like Hellschreiber and WSPR.

<div class="center border medium">

![](qrp-labs-kit-photo.jpg)

</div>

## Radio Frequency Spectral Phenomena

**Atmospheric phenomena and other special conditions can often be spotted in QRSS spectrograms**. One of the most common special cases are radio frequency reflections off of airplanes resulting in the radio waves arriving at the receiver simultaneously via two different paths (a form of [multipath propagation](https://en.wikipedia.org/wiki/Multipath_propagation)). Due to the Doppler shift from the airplane approaching the receiver the signal from the reflected path appears higher frequency than the direct path, and as the airplane flies over and begins heading away the signal from the reflected path decreases in frequency relative to the signal of the direct path. The image below is one of my favorites, captured by [Andy (G0FTD)](https://sites.google.com/site/andyg0ftd) in the 10m QRSS band. [QRSS de W4HBK](http://pensacolasnapper.blogspot.com/) is a website that has many blog posts about rare and special grabs, demonstrating effects of meteors and coronal mass ejections on QRSS signals.

<div class="center medium">

![](rf-reflection-airplane.png)

</div>

<div class="center medium border">

![](rf-reflection-airplane.jpg)

</div>

## QRSS Transmitters are Not Beacons

Radio beacons send continuous, automated, unattended, one-way transmissions without specific reception targets. In contrast, QRSS transmitters are only intended to be transmitting when the control operator is available to control them, and the recipients are known QRSS grabbers around the world. To highlight the distinction from radio beacons, QRSS transmitters are termed ***Manned Experimental Propagation Transmitters (MEPTs)***. Users in the United States will recall that the FCC (in Part 97.203) confines operation of radio beacons to specific regions of the radio spectrum and disallows operation of beacons below 28 MHz. Note that amateur radio beacons typically operate up to 100 W which is a power level multiple orders of magnitude greater than QRSS transmitters. MEPTs, in contrast, can transmit in any portion of the radio frequency spectrum where CW operation is permitted.

## The New Age of QRSS
**QRSS was first mentioned in epsisode 28 of [The Soldersmoke Podcast](http://www.soldersmoke.com/) on July 30, 2006.** It was discussed in several episodes over the next few years, and a 2009 [post about QRSS](https://hackaday.com/2009/02/22/qrss-radio-amateurs-slow-speed-narrowband/) on Hack-A-Day brought it to my attention. In the early days of QRSS the only way to transmit QRSS was to design and build your own transmitter. [David Hassall (WA5DJJ)](http://www.zianet.com/dhassall/QRSS_A.html), [Bill Houghton (W4HBK)](http://pensacolasnapper.blogspot.com/), [Hans Summers (G0UPL)](http://www.hanssummers.com/), and others would post their designs on their personal websites along with notes about where their transmitters had been spotted. In the following years the act of creating QRSS grabbers became streamlined, and websites like I2NDT's [QRSS Grabber Compendium](https://digilander.libero.it/i2ndt/grabber/grabber-compendium.htm) and [QRSS Plus](https://swharden.com/qrss/plus/) made it easier to see QRSS signals around the world. [Hans Summers](http://www.hanssummers.com/) (G0UPL) began selling QRSS transmitter kits at amateur radio conventions, then later through the [QRP Labs website](https://www.qrp-labs.com/). As more people started selling and buying kits (and documenting their experiences) it became easier and easier to get started with QRSS. Before QRSS kits were easy to obtain the only way to participate in the hobby was to design and build a transmitter from scratch, representing a high barrier to entry for those potentially interested in this fascinating hobby. Now with the availability of high quality QRSS transmitter kits and the ubiquity of internet tools and software to facilitate QRSS reception, it's easier than ever to get involved in this exciting field! For these reasons I believe we have entered into a _New Age of QRSS_.

## QRSS Frequency Bands

**This table shows the QRSS frequency range for every major amateur radio band.** Primary QRSS band windows are 100 Hz wide, although an extra 100 Hz above and below the band is typically included in grabber spectrograms to aid experimentation in these regions. [WSPR](http://wsprnet.org/) bands are often just above the QRSS window, so WSPR transmissions frequently appear on QRSS grabs.

<div class="center">

Band | Lower Edge (Hz) | Upper Edge (Hz) | Dial Frequency (Hz)
---|---|---|---
160M | 1,843,000 | 1,843,100 | 1,836,600
80M | 3,500,800 | 3,500,900 | 3,568,600
40M | 7,000,800 | 7,000,900 | 7,038,600
30M | 10,140,000 | 10,140,100 | 10,138,700
20M | 14,000,800 | 14,000,900 | 14,095,600
17M | 18,068,800 | 18,068,900 | 18,104,600
15M | 21,000,800 | 21,000,900 | 21,094,600
12M | 24,860,800 | 24,860,900 | 24,924,600
10M | 28,332,800 | 28,332,900 | 28,124,600

</div>


> **⚠️ WARNING:** These frequencies sometimes change based upon community discussion. The latest frequencies can be found at https://groups.io/g/qrssknights/wiki/3964

**When tuning your radio your dial frequency may be lower than the QRSS frequency.** If you are using upper-sideband (USB) mode, you have to tune your radio dial 1.4 kHz _below_ the QRSS band to hear QRSS signals as a 1.4 kHz tone. Recommended dial frequencies in the table above are suitable for receiving QRSS and WSPR.

## QRSS Knights

**The [QRSS Knights](https://groups.io/g/qrssknights) is a group of QRSS enthusiasts** who coordinate events and discuss experiments over email. The group is kind and welcoming to newcomers, and those interested in learning more about QRSS are encouraged to join the mailing list.

## Resources

* This page can be found at http://swharden.com/qrss

* [QRSS Plus](https://swharden.com/qrss/plus/) is an automatically-updating active QRSS grabber list

* [What is QRSS](https://www.qsl.net/m0ayf/What-is-QRSS.html) by M0AYF is a classic summary of QRSS, with many links to detailed schematics and design consideration notes. Notably the sections for [receiving QRSS](https://www.qsl.net/m0ayf/Receiving-QRSS.html) and [transmitting QRSS](https://www.qsl.net/m0ayf/Transmitting-QRSS.html) are great places to learn more.

* [QRSS and You](http://www.ka7oei.com/qrss1.html) by KA7OEI is another classic summary of QRSS.

* Weak Signal Propagation Reporter (WSPR) is a low power radio protocol that typically operates adjacent to the QRSS bands and provides automated decoding of callsign, power, and location information. Read more at http://wsprnet.org

* [The QRSS Adventure](http://www.zianet.com/dhassall/QRSS_A.html) by Dave Hassall (WA5DJJ) has circuit designs and commentary spanning far back into the early days of QRSS. His [ 1,164,000,000 Miles per Watt Test](http://www.zianet.com/dhassall/BILLIONMPW.html) is extraordinary!

* [QRSS de W4HBK](http://pensacolasnapper.blogspot.com/) website by Bill Houghton (W4HBK) contains many useful blog posts about advanced QRSS topics. The website also has many examples of special grabs depicting rare events and atmospheric phenomena.

* [Hans Summers' website](http://www.hanssummers.com/) (the founder of [QRP Labs](https://www.qrp-labs.com/)) has many excellent resources related to RF design and early work in the QRSS space.

* [Simple QRP Equipment](https://www.qsl.net/pa2ohh/) by Onno  (PA2OHH) is a collection of fantastic resources related to QRSS transmission, reception, and software design.

* [Electronics & HAM Radio Blog](http://wa0uwh.blogspot.com/) by Eldon Brown (WA0UWH) has many fantastic articles about QRSS. Eldon's SMT band-edge transmitter inspired me to make a SMT QRSS transmitter many years later.

* [Dave Richards, AA7EE](https://aa7ee.wordpress.com/) has a fantastic website documenting many amateur radio topics including QRSS. This website has the prettiest pictures of circuit boards you'll ever see.

* My [QRSS Hardware](https://github.com/swharden/QRSS-hardware) GitHub page collects notes and resources related to QRSS transmitter and receiver design.

<div class="center border">

![](smt-qrss.jpg)

</div>
Source code for ECG Simulator Circuit
---
title: ECG Simulator Circuit
date: 2020-09-27 17:11:00
tags:
  - circuit
  - diyECG
---

# ECG Simulator Circuit

**This page describes a simple circuit which produces ECG-like waveform.** The waveform is not very detailed, but it contains a sharp depolarizing (rising) component, a slower hyperpolarizing (falling) component, and a repetition rate of approximately one beat per second making it potentially useful for testing heartbeat detection circuitry.

**In 2019 I released a [YouTube video](https://www.youtube.com/watch?v=sP_-f5nsOEo) and [blog post](https://swharden.com/blog/2019-03-15-sound-card-ecg-with-ad8232/) showing how to build an ECG machine** using an AD8232 interfaced to a computer's sound card. At the end of the video I discussed how to use a 555 timer to create a waveform roughly like an ECG signal, but I didn't post the circuit at the end of that video. I get questions about it from time to time, so I'll share my best guess at what that circuit was here using LTSpice to simulate it.

<div class="center border">

![](ltspice-ecg-simulator.png)

</div>

## Design Notes

* The 555 timer generates pulses about once per second.

* The diode (D1) causes the 555 to produce very short pulses. The duty of the pulses is controlled by the resistance in series with the diode (R3), with higher resistances resulting in larger duty.

* The main purpose of the first op-amp is to invert polarity of the signal emitted by the 555. The signal is a square wave at about 1Hz, but it is mostly high with brief low pulses.

* The second op-amp serves as a voltage buffer to stabilize the output, and the final series capacitor shifts the voltage so it's centered around zero.

* Unity gain op-amps should have some feedback resistance to improve small-signal stability in production applications, but for messing around here I felt fine omitting them.

## Resources

* LTSpice file for this project: [ecg.asc](ecg.asc)

* You will need the LM741 model found on the [Using MOD Files in LTSpice](https://swharden.com/blog/2020-09-26-ltspice-mod-files/) page

* My [Action Potential Generator Circuit](https://swharden.com/blog/2017-08-12-analog-action-potential-generator-circuit/) and [Microcontroller Action Potential Generator](https://swharden.com/blog/2017-08-20-microcontroller-action-potential-generator/) articles describe method to produce a similar waveform (designed to look more like what firing neurons produce) using transistors to charge/discharge a capacitor rather than op-amps.
Source code for Using MOD Files in LTSpice
---
title: Using MOD Files in LTSpice
date: 2020-09-27 16:21:00
tags:
  - circuit
---

# Using MOD Files in LTSpice

**This page shows how to use the LM741 op-amp model file in LTSpice.** This is surprisingly un-intuitive, but is a good thing to know how to do. Model files can often be downloaded by vendor sites, but LTSpice only comes pre-loaded with models of common LT components.

## Step 1: Download a Model (.mod) File

I found [`LM741.MOD`](LM741.MOD) available on the TI's [LM741 product page](https://www.ti.com/product/LM741).

Save it wherever you want, but you will need to know the full path to this file later.

## Step 2: Determine the Name

Open the model file in a text editor and look for the line starting with `.SUBCKT`. The top of LM741.MOD looks like this:

```c
* connections:      non-inverting input
*                   |   inverting input
*                   |   |   positive power supply
*                   |   |   |   negative power supply
*                   |   |   |   |   output
*                   |   |   |   |   |
*                   |   |   |   |   |
.SUBCKT LM741/NS    1   2  99  50  28
```

The last line tells us the name of this model's sub-circuit is `LM741/NS`

## Step 3: Include the Model File

Click the ".op" button on the toolbar, then add `.include` followed by the full path to the model file. After clicking OK place the text somewhere on your LTSpice circuit diagram.

<div class="center border">

![](op2.png)

</div>

## Step 4: Insert a General Purpose Part

We know the part we are including is a 5-pin op-amp, so we can start by placing a generic component. Notice the description says _you must give the value a name and include this file_. We will do this in the next step.

<div class="center border">

![](opamp.png)

</div>

## Step 5: Configure the Component to use the Model

Right-click the op-amp and update its `Value` to match the name of the subcircuit we read from the model file earlier.

<div class="center border">

![](name.png)

</div>

## Step 6: Simulate Your Circuit

Your new component will run using the properties of the model you downloaded.

<div class="center border">

![](ltspice-lm741.png)

</div>
Source code for Exponential Fit with Python
---
title: Exponential Fit with Python
date: 2020-09-24 17:45:00
tags:
  - python
---

# Exponential Fit with Python

**Fitting an exponential curve to data is a common task** and in this example we'll use Python and SciPy to determine parameters for a curve fitted to arbitrary X/Y points. You can follow along using the [fit.ipynb](fit.ipynb) Jupyter notebook.

```python
import numpy as np
import scipy.optimize
import matplotlib.pyplot as plt

xs = np.arange(12) + 7
ys = np.array([304.08994, 229.13878, 173.71886, 135.75499,
               111.096794, 94.25109, 81.55578, 71.30187, 
               62.146603, 54.212032, 49.20715, 46.765743])

plt.plot(xs, ys, '.')
plt.title("Original Data")
```

<div class="center">

![](original.png)

</div>

**To fit an arbitrary curve** we must first define it as a function. We can then call `scipy.optimize.curve_fit` which will tweak the arguments (using arguments we provide as the starting parameters) to best fit the data. In this example we will use a single [exponential decay](https://en.wikipedia.org/wiki/Exponential_decay) function. 

```python
def monoExp(x, m, t, b):
    return m * np.exp(-t * x) + b
```

**In biology / electrophysiology _biexponential_ functions are often used** to separate fast and slow components of exponential decay which may be caused by different mechanisms and occur at different rates. In this example we will only fit the data to a method with a exponential component (a _monoexponential_ function), but the idea is the same.

```python
# perform the fit
p0 = (2000, .1, 50) # start with values near those we expect
params, cv = scipy.optimize.curve_fit(monoExp, xs, ys, p0)
m, t, b = params
sampleRate = 20_000 # Hz
tauSec = (1 / t) / sampleRate

# plot the results
plt.plot(xs, ys, '.', label="data")
plt.plot(xs, monoExp(xs, m, t, b), '--', label="fitted")
plt.title("Fitted Exponential Curve")

# inspect the parameters
print(f"Y = {m} * e^(-{t} * x) + {b}")
print(f"Tau = {tauSec * 1e6} µs")
```

<div class="center">

![](fitted.png)

</div>

```
Y = 2666.499 * e^(-0.332 * x) + 42.494
Tau = 150.422 µs
```

## Extrapolating the Fitted Curve

**We can use the calculated parameters to extend this curve** to any position by passing X values of interest into the function we used during the fit. 

**The value at time 0** is simply `m + b` because the exponential component becomes e^(0) which is 1.

```python
xs2 = np.arange(25)
ys2 = monoExp(xs2, m, t, b)

plt.plot(xs, ys, '.', label="data")
plt.plot(xs2, ys2, '--', label="fitted")
plt.title("Extrapolated Exponential Curve")
```

<div class="center">

![](fitted2.png)

</div>

## Constraining the Infinite Decay Value

**What if we know our data decays to 0?** It's not best to fit to an exponential decay function that lets the `b` component be whatever it wants. Indeed, our fit from earlier calculated the ideal `b` to be `42.494` but what if we know it should be `0`? The solution is to fit using an exponential function where `b` is constrained to 0 (or whatever value you know it to be).

```python
def monoExpZeroB(x, m, t):
    return m * np.exp(-t * x)

# perform the fit using the function where B is 0
p0 = (2000, .1) # start with values near those we expect
paramsB, cv = scipy.optimize.curve_fit(monoExpZeroB, xs, ys, p0)
mB, tB = paramsB
sampleRate = 20_000 # Hz
tauSec = (1 / tB) / sampleRate

# inspect the results
print(f"Y = {mB} * e^(-{tB} * x)")
print(f"Tau = {tauSec * 1e6} µs")

# compare this curve to the original
ys2B = monoExpZeroB(xs2, mB, tB)
plt.plot(xs, ys, '.', label="data")
plt.plot(xs2, ys2, '--', label="fitted")
plt.plot(xs2, ys2B, '--', label="zero B")
```

```
Y = 1245.580 * e^(-0.210 * x)
Tau = 237.711 µs
```

<div class="center">

![](fits.png)

</div>

**The curves produced are very different** at the extremes (especially when time is 0), even though they appear to both fit the data points nicely. Which curve is more accurate? That depends on your application. A hint can be gained by inspecting the time constants of these two curves.

<div class="center">

Parameter | Fitted B | Fixed B
---|---|---
m|2666.499|1245.580
t|0.332|0.210
Tau|150.422 µs|237.711 µs
b|42.494|0

</div>

**By inspecting Tau** I can gain insight into which method may be better for me to use in my application. I expect Tau to be near 250 µs, leading me to trust the fixed-B method over the fitted B method. Choosing the correct method has great implications on the value of `m` (which is also the value of the curve when time is 0).
page 1, page 2, page 3, page 4, page 5, page 6, page 7, page 8, page 9, page 10, page 11, page 12, page 13, page 14, page 15, page 16, page 17, page 18, page 19, page 20, page 21, page 22, page 23, page 24, page 25, page 26, page 27, page 28, page 29, page 30, page 31, page 32, page 33, page 34, page 35, page 36, page 37, page 38, page 39, page 40, page 41, page 42, page 43, page 44, page 45, page 46, page 47, page 48, page 49
All Blog Posts