The personal website of Scott W Harden

Draw Animated Graphics in the Browser with Blazor WebAssembly

Client-side Blazor allows C# graphics models to be rendered in the browser. This means .NET developers can create web apps with business logic written (and tested) in C# instead of being forced to write their business logic in a totally different language (JavaScript) just to get it to run in the browser. In this project we'll create an interactive graphics model (a field of balls that bounce off the edge of the screen) entirely in C#, and draw the model on the screen using an API that allows us to interact with a HTML5 Canvas.

Strategy

At the time of this writing Blazor WebAssembly can't directly paint on the screen, so JavaScript is required somewhere to make this happen. The Blazor.Extensions.Canvas package has a Canvas component and provides a C# API for all of its JavaScript methods, allowing you to draw on the canvas entirely from C#.

Rendering gets a little slower every time a JavaScript function is called from Blazor. If a large number of shapes are required (a lot of JavaScript calls), performance may be limited compared to a similar application written entirely in JavaScript. For high performance rendering of many objects, a rendering loop written entirely in JavaScript may be required.

This method is good for simple models with a limited number of shapes. It allows the business logic (and tests) to remain entirely in C#. The same graphics model code could be displayed in the browser (using HTML canvas) or on the desktop in WPF and WinForms apps (using System.Drawing or SkiaSharp).

Step 1: Create a Pure C# Graphics Model

I'm using graphics model to describe the state of the field of balls and logic required to move each ball around. Ideally this logic would be isolated in a separate library. At the time of writing a .NET Standard C# library seems like a good idea, so the same graphics model could be used in .NET Framework and .NET Core environments.

Models/Ball.cs

public class Ball
{
    public double X { get; private set; }
    public double Y { get; private set; }
    public double XVel { get; private set; }
    public double YVel { get; private set; }
    public double Radius { get; private set; }
    public string Color { get; private set; }

    public Ball(double x, double y, double xVel, double yVel, double radius, string color)
    {
        (X, Y, XVel, YVel, R, Color) = (x, y, xVel, yVel, radius, color);
    }

    public void StepForward(double width, double height)
    {
        X += XVel;
        Y += YVel;
        if (X < 0 || X > width)
            XVel *= -1;
        if (Y < 0 || Y > height)
            YVel *= -1;

        if (X < 0)
            X += 0 - X;
        else if (X > width)
            X -= X - width;

        if (Y < 0)
            Y += 0 - Y;
        if (Y > height)
            Y -= Y - height;
    }
}

Models/Field.cs

public class Field
{
    public readonly List<Ball> Balls = new List<Ball>();
    public double Width { get; private set; }
    public double Height { get; private set; }

    public void Resize(double width, double height) =>
        (Width, Height) = (width, height);

    public void StepForward()
    {
        foreach (Ball ball in Balls)
            ball.StepForward(Width, Height);
    }

    private double RandomVelocity(Random rand, double min, double max)
    {
        double v = min + (max - min) * rand.NextDouble();
        if (rand.NextDouble() > .5)
            v *= -1;
        return v;
    }


    private string RandomColor(Random rand) => 
        string.Format("#{0:X6}", rand.Next(0xFFFFFF));

    public void AddRandomBalls(int count = 10)
    {
        double minSpeed = .5;
        double maxSpeed = 5;
        double radius = 10;
        Random rand = new Random();

        for (int i = 0; i < count; i++)
        {
            Balls.Add(
                new Ball(
                    x: Width * rand.NextDouble(),
                    y: Height * rand.NextDouble(),
                    xVel: RandomVelocity(rand, minSpeed, maxSpeed),
                    yVel: RandomVelocity(rand, minSpeed, maxSpeed),
                    radius: radius,
                    color: RandomColor(rand);
                )
            );
        }
    }
}

Step 2: Get the Blazor.Extensions.Canvas Package

Use NuGet to install Blazor.Extensions.Canvas

Install-Package Blazor.Extensions.Canvas

Step 3: Add a script to index.html

This JavaScript sets-up the render loop which automatically calls RenderInBlazor C# method on every frame. It also calls the ResizeInBlazor C# function whenever the canvas is resized so the graphics model's dimensions can be updated. You can place it just before the closing </body> tag.

This code will automatically resize the canvas and graphics model to fit the window, but for a fixed-size canvas you can omit the addEventListener and resizeCanvasToFitWindow calls.

<script src='_content/Blazor.Extensions.Canvas/blazor.extensions.canvas.js'></script>
<script>
    function renderJS(timeStamp) {
        theInstance.invokeMethodAsync('RenderInBlazor', timeStamp);
        window.requestAnimationFrame(renderJS);
    }

    function resizeCanvasToFitWindow() {
        var holder = document.getElementById('canvasHolder');
        var canvas = holder.querySelector('canvas');
        if (canvas) {
            canvas.width = window.innerWidth;
            canvas.height = window.innerHeight;
            theInstance.invokeMethodAsync('ResizeInBlazor', canvas.width, canvas.height);
        }
    }

    window.initRenderJS = (instance) => {
        window.theInstance = instance;
        window.addEventListener("resize", resizeCanvasToFitWindow);
        resizeCanvasToFitWindow();
        window.requestAnimationFrame(renderJS);
    };
</script>

Step 4: Create a page for the Canvas and Render code

This step ties everything together:

  • The graphics model: a private class stored at this level
  • The canvas component: a protected field
  • The canvas: a Razor component referencing the canvas component
  • The init code: which tells JavaScipt to start the render loop
  • The render method: C# function called from the JavaScript render loop when a frame is to be drawn
  • The resize method: C# function called from JavaScript to update the model when the canvas size changes
  • JavaScript runtime injection: This allows Blazor/JavaScript interoperability (JS interop)

For simplicity it's demonstrated here using a code-behind, but a clearer strategy would be to move the render logic into its own class/file.

@page "/"

@using Blazor.Extensions
@using Blazor.Extensions.Canvas
@using Blazor.Extensions.Canvas.Canvas2D
@inject IJSRuntime JsRuntime;

<div id="canvasHolder" style="position: fixed; width: 100%; height: 100%">
    <BECanvas Width="600" Height="400" @ref="CanvasRef"></BECanvas>
</div>

@code{
    private Models.Field BallField = new Models.Field();
    private Canvas2DContext ctx;
    protected BECanvasComponent CanvasRef;
    private DateTime LastRender;

    protected override async Task OnAfterRenderAsync(bool firstRender)
    {
        this.ctx = await CanvasRef.CreateCanvas2DAsync();
        await JsRuntime.InvokeAsync<object>("initRenderJS", DotNetObjectReference.Create(this));
        await base.OnInitializedAsync();
    }

    [JSInvokable]
    public void ResizeInBlazor(double width, double height) => BallField.Resize(width, height);

    [JSInvokable]
    public async ValueTask RenderInBlazor(float timeStamp)
    {
        if (BallField.Balls.Count == 0)
            BallField.AddRandomBalls(50);
        BallField.StepForward();

        double fps = 1.0 / (DateTime.Now - LastRender).TotalSeconds;
        LastRender = DateTime.Now;

        await this.ctx.BeginBatchAsync();
        await this.ctx.ClearRectAsync(0, 0, BallField.Width, BallField.Height);
        await this.ctx.SetFillStyleAsync("#003366");
        await this.ctx.FillRectAsync(0, 0, BallField.Width, BallField.Height);
        await this.ctx.SetFontAsync("26px Segoe UI");
        await this.ctx.SetFillStyleAsync("#FFFFFF");
        await this.ctx.FillTextAsync("Blazor WebAssembly + HTML Canvas", 10, 30);
        await this.ctx.SetFontAsync("16px consolas");
        await this.ctx.FillTextAsync($"FPS: {fps:0.000}", 10, 50);
        await this.ctx.SetStrokeStyleAsync("#FFFFFF");
        foreach (var ball in BallField.Balls)
        {
            await this.ctx.BeginPathAsync();
            await this.ctx.ArcAsync(ball.X, ball.Y, ball.Radius, 0, 2 * Math.PI, false);
            await this.ctx.SetFillStyleAsync(ball.Color);
            await this.ctx.FillAsync();
            await this.ctx.StrokeAsync();
        }
        await this.ctx.EndBatchAsync();
    }
}

Notice how the JavaScript calls are wrapped in BeginBatchAsync and EndBatchAsync. This allows the calls between these two statements to be bundled into a single call. Since Blazor/JavaScript interop calls are the primary bottleneck in this system, limiting the number of individual JavaScript calls has a large influence on final performance.

Notes and References

💡 I observed a strong performance increase when upgrading from .NET Core 3.1 to .NET 5. If you can, Migrate your Blazor project to .NET 5 or newer

Markdown source code last modified on January 18th, 2021
---
title: Draw Animated Graphics in the Browser with Blazor WebAssembly
date: 2021-01-07 21:46:00
tags: blazor, csharp
---

# Draw Animated Graphics in the Browser with Blazor WebAssembly

**Client-side Blazor allows C# graphics models to be rendered in the browser.** This means .NET developers can create web apps with business logic written (and tested) in C# instead of being forced to write their business logic in a totally different language (JavaScript) just to get it to run in the browser. In this project we'll create an interactive graphics model (a field of balls that bounce off the edge of the screen) entirely in C#, and draw the model on the screen using an API that allows us to interact with a HTML5 Canvas.

<div class="text-center">

[![](blazor-canvas-demo.gif)](app)

</div>

* [**View the live demo**](app)
* [**Download the source code**](blazor-canvas.zip)

## Strategy

**At the time of this writing Blazor WebAssembly can't directly paint on the screen, so JavaScript is required somewhere to make this happen.** The [`Blazor.Extensions.Canvas`](https://github.com/BlazorExtensions/Canvas) package has a Canvas component and provides a C# API for all of its JavaScript methods, allowing you to draw on the canvas entirely from C#. 

**Rendering gets a little slower every time a JavaScript function is called from Blazor.** If a large number of shapes are required (a lot of JavaScript calls), performance may be limited compared to a similar application written entirely in JavaScript. For high performance rendering of many objects, a rendering loop written entirely in JavaScript may be required.

**This method is good for simple models with a limited number of shapes.** It allows the business logic (and tests) to remain entirely in C#. The same graphics model code could be displayed in the browser (using HTML canvas) or on the desktop in WPF and WinForms apps (using System.Drawing or SkiaSharp).

## Step 1: Create a Pure C# Graphics Model

**I'm using _graphics model_ to describe the state of the field of balls and logic required to move each ball around.** Ideally this logic would be isolated in a separate library. At the time of writing a .NET Standard C# library seems like a good idea, so the same graphics model could be used in .NET Framework and .NET Core environments.

### Models/Ball.cs

```cs
public class Ball
{
    public double X { get; private set; }
    public double Y { get; private set; }
    public double XVel { get; private set; }
    public double YVel { get; private set; }
    public double Radius { get; private set; }
    public string Color { get; private set; }

    public Ball(double x, double y, double xVel, double yVel, double radius, string color)
    {
        (X, Y, XVel, YVel, R, Color) = (x, y, xVel, yVel, radius, color);
    }

    public void StepForward(double width, double height)
    {
        X += XVel;
        Y += YVel;
        if (X < 0 || X > width)
            XVel *= -1;
        if (Y < 0 || Y > height)
            YVel *= -1;

        if (X < 0)
            X += 0 - X;
        else if (X > width)
            X -= X - width;

        if (Y < 0)
            Y += 0 - Y;
        if (Y > height)
            Y -= Y - height;
    }
}
```

### Models/Field.cs

```cs
public class Field
{
    public readonly List<Ball> Balls = new List<Ball>();
    public double Width { get; private set; }
    public double Height { get; private set; }

    public void Resize(double width, double height) =>
        (Width, Height) = (width, height);

    public void StepForward()
    {
        foreach (Ball ball in Balls)
            ball.StepForward(Width, Height);
    }

    private double RandomVelocity(Random rand, double min, double max)
    {
        double v = min + (max - min) * rand.NextDouble();
        if (rand.NextDouble() > .5)
            v *= -1;
        return v;
    }


    private string RandomColor(Random rand) => 
	    string.Format("#{0:X6}", rand.Next(0xFFFFFF));
	
    public void AddRandomBalls(int count = 10)
    {
        double minSpeed = .5;
        double maxSpeed = 5;
        double radius = 10;
        Random rand = new Random();

        for (int i = 0; i < count; i++)
        {
            Balls.Add(
                new Ball(
                    x: Width * rand.NextDouble(),
                    y: Height * rand.NextDouble(),
                    xVel: RandomVelocity(rand, minSpeed, maxSpeed),
                    yVel: RandomVelocity(rand, minSpeed, maxSpeed),
                    radius: radius,
                    color: RandomColor(rand);
                )
            );
        }
    }
}
```

## Step 2: Get the Blazor.Extensions.Canvas Package

Use NuGet to install [Blazor.Extensions.Canvas](https://github.com/BlazorExtensions/Canvas)

```cs
Install-Package Blazor.Extensions.Canvas
```

## Step 3: Add a script to index.html

**This JavaScript sets-up the render loop** which automatically calls `RenderInBlazor` C# method on every frame. It also calls the `ResizeInBlazor` C# function whenever the canvas is resized so the graphics model's dimensions can be updated. You can place it just before the closing `</body>` tag.

**This code will automatically resize the canvas and graphics model to fit the window,** but for a fixed-size canvas you can omit the `addEventListener` and `resizeCanvasToFitWindow` calls.

```html
<script src='_content/Blazor.Extensions.Canvas/blazor.extensions.canvas.js'></script>
<script>
    function renderJS(timeStamp) {
        theInstance.invokeMethodAsync('RenderInBlazor', timeStamp);
        window.requestAnimationFrame(renderJS);
    }

    function resizeCanvasToFitWindow() {
        var holder = document.getElementById('canvasHolder');
        var canvas = holder.querySelector('canvas');
        if (canvas) {
            canvas.width = window.innerWidth;
            canvas.height = window.innerHeight;
            theInstance.invokeMethodAsync('ResizeInBlazor', canvas.width, canvas.height);
        }
    }

    window.initRenderJS = (instance) => {
        window.theInstance = instance;
        window.addEventListener("resize", resizeCanvasToFitWindow);
        resizeCanvasToFitWindow();
        window.requestAnimationFrame(renderJS);
    };
</script>
```

## Step 4: Create a page for the Canvas and Render code

This step ties everything together:
 * **The graphics model**: a private class stored at this level
 * **The canvas component**: a protected field
 * **The canvas**: a Razor component referencing the canvas component
 * **The init code**: which tells JavaScipt to start the render loop
 * **The render method**: C# function called from the JavaScript render loop when a frame is to be drawn
 * **The resize method**: C# function called from JavaScript to update the model when the canvas size changes
* **JavaScript runtime injection**: This allows Blazor/JavaScript interoperability (JS interop)

For simplicity it's demonstrated here using a code-behind, but a clearer strategy would be to move the render logic into its own class/file.

```cs
@page "/"

@using Blazor.Extensions
@using Blazor.Extensions.Canvas
@using Blazor.Extensions.Canvas.Canvas2D
@inject IJSRuntime JsRuntime;

<div id="canvasHolder" style="position: fixed; width: 100%; height: 100%">
    <BECanvas Width="600" Height="400" @ref="CanvasRef"></BECanvas>
</div>

@code{
    private Models.Field BallField = new Models.Field();
    private Canvas2DContext ctx;
    protected BECanvasComponent CanvasRef;
    private DateTime LastRender;

    protected override async Task OnAfterRenderAsync(bool firstRender)
    {
        this.ctx = await CanvasRef.CreateCanvas2DAsync();
        await JsRuntime.InvokeAsync<object>("initRenderJS", DotNetObjectReference.Create(this));
        await base.OnInitializedAsync();
    }

    [JSInvokable]
    public void ResizeInBlazor(double width, double height) => BallField.Resize(width, height);

    [JSInvokable]
    public async ValueTask RenderInBlazor(float timeStamp)
    {
        if (BallField.Balls.Count == 0)
            BallField.AddRandomBalls(50);
        BallField.StepForward();

        double fps = 1.0 / (DateTime.Now - LastRender).TotalSeconds;
        LastRender = DateTime.Now;

        await this.ctx.BeginBatchAsync();
        await this.ctx.ClearRectAsync(0, 0, BallField.Width, BallField.Height);
        await this.ctx.SetFillStyleAsync("#003366");
        await this.ctx.FillRectAsync(0, 0, BallField.Width, BallField.Height);
        await this.ctx.SetFontAsync("26px Segoe UI");
        await this.ctx.SetFillStyleAsync("#FFFFFF");
        await this.ctx.FillTextAsync("Blazor WebAssembly + HTML Canvas", 10, 30);
        await this.ctx.SetFontAsync("16px consolas");
        await this.ctx.FillTextAsync($"FPS: {fps:0.000}", 10, 50);
        await this.ctx.SetStrokeStyleAsync("#FFFFFF");
        foreach (var ball in BallField.Balls)
        {
            await this.ctx.BeginPathAsync();
            await this.ctx.ArcAsync(ball.X, ball.Y, ball.Radius, 0, 2 * Math.PI, false);
            await this.ctx.SetFillStyleAsync(ball.Color);
            await this.ctx.FillAsync();
            await this.ctx.StrokeAsync();
        }
        await this.ctx.EndBatchAsync();
    }
}
```

**Notice how the JavaScript calls are wrapped in `BeginBatchAsync` and `EndBatchAsync`.** This allows the calls between these two statements to be bundled into a single call. Since Blazor/JavaScript interop calls are the primary bottleneck in this system, limiting the number of individual JavaScript calls has a large influence on final performance.

## Notes and References

> **💡 I observed a strong performance increase when upgrading from .NET Core 3.1 to .NET 5.** If you can, [Migrate your Blazor project to .NET 5](https://docs.microsoft.com/en-us/aspnet/core/migration/31-to-50#update-blazor-webassembly-projects) or newer

* **A hybrid JavaScript/C# architecture is be possible** involving exchange of graphics model _data_, whereby JavaScript is used to render a model but C# is used to advanced the model. This could be achieved with a single interop call passing model data as JSON. This is explored in the later blog posts [Mystify your Browser with Blazor](https://swharden.com/blog/2021-01-09-blazor-mystify/) and [Boids in your Browser with Blazor](https://swharden.com/blog/2021-01-08-blazor-boids/).

* **Live demo of this project**: [Ball field demo](app)

* **Source code for this project**: [blazor-canvas.zip](app)

* **C# Data Visualization:** https://github.com/swharden/Csharp-Data-Visualization

Create a new .NET 5 Blazor App with Visual Studio 2019

When Visual Studio 2019 creates a new Blazor app it defaults to .NET Standard 2.1. A few months ago .NET 5 was released, and upgrading Blazor apps to use .NET 5 has many strong performance improvements. I'm sure a Visual Studio update will soon make this easier, but for now to create a new .NET Blazor app running .NET 5 I use these steps:

  • Create a new .NET core 3.1 Blazor App
  • Edit the csproj file
    • Change SDK from Microsoft.NET.Sdk.Web to Microsoft.NET.Sdk.BlazorWebAssembly
    • Remove RazorLangVersion
    • Update TargetFramework from netstandard2.1 to net5.0
    • Remove the Microsoft.AspNetCore.Components.WebAssembly.Build package reference
  • update all NuGet packages to their latest versions

Extensive details can be found on Microsoft's official Migrate from ASP.NET Core 3.1 to 5.0 documentation page, but I find this short list of steps easier to refer to.

Download a New .NET 5.0.1 Blazor App

Here is a new project running .NET 5.0.1 with Bootstrap 5.0 alpha:

To change the project name/namespace:

  • Rename the .sln and .csproj files
  • Update the .sln file in a text editor to match the new filenames
  • Update the namespace in Program.cs
  • Update the namespaces in _Imports.razor
Markdown source code last modified on January 18th, 2021
---
title: Create a new .NET 5 Blazor Project
date: 2021-01-06 19:39:00
tags: blazor, csharp
---

## Create a new .NET 5 Blazor App with Visual Studio 2019

**When Visual Studio 2019 creates a new Blazor app it defaults to .NET Standard 2.1**. A few months ago .NET 5 was released, and upgrading Blazor apps to use .NET 5 has many strong performance improvements. I'm sure a Visual Studio update will soon make this easier, but for now to create a new .NET Blazor app running .NET 5 I use these steps:

* Create a new .NET core 3.1 Blazor App
* Edit the csproj file
  * Change SDK from `Microsoft.NET.Sdk.Web` to `Microsoft.NET.Sdk.BlazorWebAssembly`
  * Remove `RazorLangVersion` 
  * Update `TargetFramework` from `netstandard2.1` to `net5.0`
  * Remove the `Microsoft.AspNetCore.Components.WebAssembly.Build` package reference
* update all NuGet packages to their latest versions

Extensive details can be found on Microsoft's official [Migrate from ASP.NET Core 3.1 to 5.0](https://docs.microsoft.com/en-us/aspnet/core/migration/31-to-50?view=aspnetcore-5.0&tabs=visual-studio#update-blazor-webassembly-projects) documentation page, but I find this short list of steps easier to refer to.

## Download a New .NET 5.0.1 Blazor App

Here is a new project running .NET 5.0.1 with Bootstrap 5.0 alpha:

* [**NewBlazorApp-net5.0.1.zip**](NewBlazorApp-net5.0.1.zip)

To change the project name/namespace:
* Rename the `.sln` and `.csproj` files
* Update the `.sln` file in a text editor to match the new filenames
* Update the namespace in `Program.cs`
* Update the namespaces in `_Imports.razor`

Show Build Date in Blazor Apps

I find it useful to add build date and .NET version to the bottom of my client-side applications. Something like:

MyApp version 1.2.3 | Built on December 29, 2020 | Running on .NET 5.0.1

I'm documenting how I do this so I can refer to it later, and also so it may be helpful to others.

@page "/"

<h1>New Blazor App</h1>
<div>App version @AppVersion</div>
<div>Running on .NET @Environment.Version</div>

@code{
    private string AppVersion
    {
        get
        {
            Version version = System.Reflection.Assembly.GetExecutingAssembly().GetName().Version;
            return $"{version.Major}.{version.Minor}.{version.Build}";
        }
    }
}

Build Date

You can't access build date entirely from code, but you can create a file containing the build date on every build then consume that file as a resource.

Step1: Add a pre-build instruction

  • Right-click the project and select "Properties"
  • Navigate to the "Build Events" section
  • Add a pre-build command: echo %date% %time% > "$(ProjectDir)\Resources\BuildDate.txt"
  • Rebuild the application

⚠️ This command assumes you are building on Windows

Step2: Add the date file as a resource

  • Right-click the project and select "Properties"
  • Navigate to the "Resources" section
  • Click "Add Resource", select "Add existing file", and choose the new text file

Step3: Reference the build date resource in code

private string BuildDateString => 
    DateTime.Parse(Properties.Resources.BuildDate).ToString("MMMM dd, yyyy");
Markdown source code last modified on January 18th, 2021
---
title: Display Build Details in Client-Size Blazor Apps
date: 2020-12-29 19:54:00
tags: csharp, blazor
---

# Show Build Date in Blazor Apps

I find it useful to add build date and .NET version to the bottom of my client-side applications. Something like:

> MyApp version 1.2.3 | Built on December 29, 2020 | Running on .NET 5.0.1

I'm documenting how I do this so I can refer to it later, and also so it may be helpful to others.

```cs
@page "/"

<h1>New Blazor App</h1>
<div>App version @AppVersion</div>
<div>Running on .NET @Environment.Version</div>

@code{
	private string AppVersion
	{
		get
		{
			Version version = System.Reflection.Assembly.GetExecutingAssembly().GetName().Version;
			return $"{version.Major}.{version.Minor}.{version.Build}";
		}
	}
}
```

## Build Date

You can't access build date entirely from code, but you _can_ create a file containing the build date on every build then consume that file as a resource.

**Step1: Add a pre-build instruction**
* Right-click the project and select "Properties"
* Navigate to the "Build Events" section
* Add a pre-build command: `echo %date% %time% > "$(ProjectDir)\Resources\BuildDate.txt"`
* Rebuild the application

> ⚠️ This command assumes you are building on Windows

**Step2: Add the date file as a resource**
* Right-click the project and select "Properties"
* Navigate to the "Resources" section
* Click "Add Resource", select "Add existing file", and choose the new text file

**Step3: Reference the build date resource in code**

```cs
private string BuildDateString => 
    DateTime.Parse(Properties.Resources.BuildDate).ToString("MMMM dd, yyyy");
```

Seven Years of QRSS Plus

This article was written for Andy (G0FTD) for publication in the December 2020 edition of 74!, The Knights QRSS Newsletter. Join the QRSS Knights mailing list for the latest QRSS news.

The QRSS hobby owes much of its success to the extraordinary efforts of amateur radio operators who run and maintain QRSS grabber stations. QRSS grabbers are built by pairing a radio receiver with a computer running software to continuously convert the received signals into spectrogram images which are uploaded to the internet every few minutes. QRSS Plus is a website that displays these spectrograms from active QRSS grabbers around the world. This article discusses the history of QRSS Plus, the technical details that make it possible, and highlights its most active contributors in 2020.

Early Days of QRSS Grabber Websites

In the early 2010s when I started to become active in the QRSS community, one of my favorite grabber websites was I2NDT's QRSS Knights Grabber Compendium. I remember checking that website from my laptop during class several times throughout the day to see if my signal was making it out of Gainesville, Florida. I also recall many nights at the bench tweaking my transmitter and looking at all the grabs from around the world trying to spot my signal.

A common problem with QRSS grabber websites was the persistance of outdated grabber images. Whenever a grabber uploaded a new spectrogram image it replaced the previous one, but when a grabber stopped uploading new images the old one would remain. Most uploaded spectrograms have a date written in small text in the corner, but at a glance (and especially in thumbnails) it was difficult to identify which grabber images were current and which were outdated.

History of QRSS Plus

I created QRSS Plus in July 2013 to solve the problem of outdated spectrograms appearing on grabber websites. QRSS Plus works by downloading grabber images every 10 minutes and recording their MD5 hash (a way to convert an image into a unique set of letters such that when the image changes the letters change). Grabbers were marked "active" if the MD5 hash from their newest image was different than the one from the previous image. This original system was entirely operated as a PHP script which ran on the back-end of a web server triggered by a cron job to download new images and update a database every 10 minutes. The primary weakness of this method was that downloading all those images took a lot of time (they were downloaded sequentially on the server). PHP is famously single-threaded, and my web host limited how long PHP scripts could run, limiting the maximum number of supported grabbers.

The back-end of QRSS Plus was redesigned in 2016 when I changed hosting companies. The new company allowed me to execute python scripts on the server, so I was no longer limited by the constraints of PHP. I redesigned QRSS Plus to download, hash, and store images every 10 minutes. This allowed QRSS Plus to display a running history of the last several grabs for each grabber, as well as support automated image stacking (averaging the last several images together to improve visualization of weak, repetitive signals). This solution is still limited by CPU time (the number of CPU seconds per day is capped by my hosting company), but continuously operating QRSS Plus does not occupy a large portion of that time.

QRSS Plus Activity in 2020

I started logging grabber updates in September 2018, allowing me to reflect on the last few years of grabber activity. It takes a lot of effort to set-up and maintain a quality QRSS grabber, and the enthusiasm and dedication of the QRSS community is impressive and inspiring!

In 2020 our community saw 155 active grabber stations! On average there were more than 60 active stations running on any given day, and the number of active stations is visibly increasing with time.

In 2020 QRSS Plus analyzed a mean of 6,041 spectrograms per day. In total, QRSS Plus analyzed over 2.2 million spectrograms this year!

This bar graph depicts the top 50 most active grabber stations ranked according to their total unique spectrogram upload count. Using this metric grabbers that update once every 10 minutes will appear to have twice as many unique grabber images as those which upload a unique image every 20 minutes.

Many QRSS grabber operators maintain multiple stations, and I think those operators deserve special attention! This year's winner for the most active contributor goes to David Hassall (WA5DJJ) who alone is responsible for 15.26% of the total number of uploaded spectrograms in 2020 🏆

The top 25 contributors with the greatest number of unique uploads in 2020 were (in order): WA5DJJ, WD4ELG, W6REK, G3VYZ, KL7L, G4IOG, W4HBK, G3YXM, HB9FXX, PA2OHH, EA8BVP, G0MQW, SA6BSS, WD4AH, 7L1RLL, VA3ROM, VE7IGH, DL4DTL, K5MO, LA5GOA, VE3GTC, AJ4VD, K4RCG, GM4SFW, and OK1FCX.

Maintaining QRSS Plus

I want to recognize Andy (G0FTD) for doing an extraordinary job maintaining the QRSS Plus website over the last several years! Although I (AJ4VD) created and maintain the QRSS Plus back-end, once it is set-up it largely operates itself. The QRSS grabber list, on the other hand, requires frequent curation. Andy has done a fantastic job monitoring the QRSS Knights mailing list and updating the grabber list in response to updates posted there so it always contains the newest grabbers and uses the latest URLs. On behalf of everyone who enjoys using QRSS Plus, thank you for your work Andy!

The Future of QRSS Plus

Today QRSS Plus is functional, but I think its front-end could be greatly improved. It is written using vanilla JavaScript, but I think moving to a front-end framework like React makes a lot of sense. Currently PHP generates HTML containing grabber data when the page is requested, but a public JSON API would make a lot of sense and make QRSS Plus it easier to develop and test. From a UX standpoint, the front-end could benefit from a simpler design that displays well on mobile and desktop browsers. I think the usability of historical grabs could be greatly improved as well. From a back-end perspective, I'd much prefer to run the application using a service like Azure or AWS rather than lean on a traditional personal website hosting plan to manage the downloads and image processing. Automatic creation of 8-hour (stitched and compressed) grabber frames seems feasible as well. It is unlikely I will work toward these improvements soon, but if you are a front-end web developer interested in working on a public open-source project, send me an email and I'd be happy to discuss how we can improve QRSS Plus together!

QRSS is a growing hobby, and if the rise in grabbers over the last few years is an indication of what the next few years will look like, I'm excited to see where the hobby continues to go! I encourage you to consider running a grabber (and to encourage others to do the same), and to continue to thank all the grabber operators and maintainers out there who make this hobby possible.

Notes and Resources

  • Data includes Jan 1 2020 through Dec 11 2020
  • Stations with <1000 unique uploads were excluded from most analyses
  • Summary data (a table of unique images per day per station) is available: qrss-plus-2020.xlsx
  • Bar graphs, scatter plots, and line charts were created with ScottPlot
  • QRSS Plus is open source on GitHub
  • A modern introduction to QRSS: The New Age of QRSS
  • FSKview is a new QRSS and WSPR Spectrogram Viewer for Windows
Markdown source code last modified on January 18th, 2021
---
title: Seven Years of QRSS Plus
date: 2020-12-14 21:16:00
tags: qrss, amateur radio
---

# Seven Years of QRSS Plus

> This article was written for [Andy (G0FTD)](https://sites.google.com/view/andy-g0ftd/) for publication in the December 2020 edition of [74!, The Knights QRSS Newsletter](https://swharden.com/qrss/74/). Join the [QRSS Knights mailing list](https://groups.io/g/qrssknights) for the latest QRSS news.

The QRSS hobby owes much of its success to the extraordinary efforts of amateur radio operators who run and maintain QRSS grabber stations. QRSS grabbers are built by pairing a radio receiver with a computer running software to continuously convert the received signals into spectrogram images which are uploaded to the internet every few minutes. [**QRSS Plus**](https://swharden.com/qrss/plus/) is a website that displays these spectrograms from active QRSS grabbers around the world. This article discusses the history of QRSS Plus, the technical details that make it possible, and highlights its most active contributors in 2020.

### Early Days of QRSS Grabber Websites

In the early 2010s when I started to become active in the QRSS community, one of my favorite grabber websites was [I2NDT's QRSS Knights Grabber Compendium](https://digilander.libero.it/i2ndt/grabber/grabber-compendium.htm). I remember checking that website from my laptop during class several times throughout the day to see if my signal was making it out of Gainesville, Florida. I also recall many nights at the bench tweaking my transmitter and looking at all the grabs from around the world trying to spot my signal. 

A common problem with QRSS grabber websites was the persistance of outdated grabber images. Whenever a grabber uploaded a new spectrogram image it replaced the previous one, but when a grabber stopped uploading new images the old one would remain. Most uploaded spectrograms have a date written in small text in the corner, but at a glance (and especially in thumbnails) it was difficult to identify which grabber images were current and which were outdated. 

### History of QRSS Plus

I created [QRSS Plus](https://swharden.com/qrss/plus/) in July 2013 to solve the problem of outdated spectrograms appearing on grabber websites. QRSS Plus works by downloading grabber images every 10 minutes and recording their [MD5 hash](https://en.wikipedia.org/wiki/MD5) (a way to convert an image into a unique set of letters such that when the image changes the letters change). Grabbers were marked "active" if the MD5 hash from their newest image was different than the one from the previous image. This original system was entirely operated as a PHP script which ran on the back-end of a web server triggered by a cron job to download new images and update a database every 10 minutes. The primary weakness of this method was that downloading all those images took a lot of time (they were downloaded sequentially on the server). PHP is famously single-threaded, and my web host limited how long PHP scripts could run, limiting the maximum number of supported grabbers.

The back-end of QRSS Plus was redesigned in 2016 when I changed hosting companies. The new company allowed me to execute python scripts on the server, so I was no longer limited by the constraints of PHP. I redesigned QRSS Plus to download, hash, and store images every 10 minutes. This allowed QRSS Plus to display a running history of the last several grabs for each grabber, as well as support automated image stacking (averaging the last several images together to improve visualization of weak, repetitive signals). This solution is still limited by CPU time (the number of CPU seconds per day is capped by my hosting company), but continuously operating QRSS Plus does not occupy a large portion of that time.

### QRSS Plus Activity in 2020

I started logging grabber updates in September 2018, allowing me to reflect on the last few years of grabber activity. It takes a lot of effort to set-up and maintain a quality QRSS grabber, and the enthusiasm and dedication of the QRSS community is impressive and inspiring! 

<div class="text-center">

![](grabbers-per-day.png)

</div>

In 2020 our community saw ***155 active grabber stations***! On average there were more than 60 active stations running on any given day, and the number of active stations is visibly increasing with time.

<div class="text-center">

![](grabs-per-day.png)

</div>

In 2020 QRSS Plus analyzed a mean of 6,041 spectrograms per day. In total, QRSS Plus analyzed over ***2.2 million spectrograms*** this year!

<div class="text-center">

![](grabbers-leader.png)

</div>

This bar graph depicts the top 50 most active grabber stations ranked according to their total unique spectrogram upload count. Using this metric grabbers that update once every 10 minutes will appear to have twice as many unique grabber images as those which upload a unique image every 20 minutes.

<div class="text-center">

![](qrss-2020-pie.png)

</div>

Many QRSS grabber operators maintain multiple stations, and I think those operators deserve special attention! This year's winner for the most active contributor goes to David Hassall (WA5DJJ) who alone is responsible for 15.26% of the total number of uploaded spectrograms in 2020 🏆

The top 25 contributors with the greatest number of unique uploads in 2020 were (in order): WA5DJJ, WD4ELG, W6REK, G3VYZ, KL7L, G4IOG, W4HBK, G3YXM, HB9FXX, PA2OHH, EA8BVP, G0MQW, SA6BSS, WD4AH, 7L1RLL, VA3ROM, VE7IGH, DL4DTL, K5MO, LA5GOA, VE3GTC, AJ4VD, K4RCG, GM4SFW, and OK1FCX. 

### Maintaining QRSS Plus

I want to recognize Andy (G0FTD) for doing an extraordinary job maintaining the QRSS Plus website over the last several years! Although I (AJ4VD) created and maintain the QRSS Plus back-end, once it is set-up it largely operates itself. The QRSS grabber list, on the other hand, requires frequent curation. Andy has done a fantastic job monitoring the [QRSS Knights mailing list](https://groups.io/g/qrssknights) and updating the grabber list in response to updates posted there so it always contains the newest grabbers and uses the latest URLs. On behalf of everyone who enjoys using QRSS Plus, thank you for your work Andy!

### The Future of QRSS Plus

Today QRSS Plus is functional, but I think its front-end could be greatly improved. It is written using vanilla JavaScript, but I think moving to a front-end framework like React makes a lot of sense. Currently PHP generates HTML containing grabber data when the page is requested, but a public JSON API would make a lot of sense and make QRSS Plus it easier to develop and test. From a UX standpoint, the front-end could benefit from a simpler design that displays well on mobile and desktop browsers. I think the usability of historical grabs could be greatly improved as well. From a back-end perspective, I'd much prefer to run the application using a service like Azure or AWS rather than lean on a traditional personal website hosting plan to manage the downloads and image processing. Automatic creation of 8-hour (stitched and compressed) grabber frames seems feasible as well. It is unlikely I will work toward these improvements soon, but if you are a front-end web developer interested in working on a public open-source project, send me an email and I'd be happy to discuss how we can improve QRSS Plus together!

QRSS is a growing hobby, and if the rise in grabbers over the last few years is an indication of what the next few years will look like, I'm excited to see where the hobby continues to go! I encourage you to consider running a grabber (and to encourage others to do the same), and to continue to thank all the grabber operators and maintainers out there who make this hobby possible.

### Notes and Resources
* Data includes Jan 1 2020 through Dec 11 2020
* Stations with <1000 unique uploads were excluded from most analyses
* Summary data (a table of unique images per day per station) is available: [qrss-plus-2020.xlsx](qrss-plus-2020.xlsx)
* Bar graphs, scatter plots, and line charts were created with [ScottPlot](https://swharden.com/scottplot)
* QRSS Plus is [open source on GitHub](https://github.com/swharden/qrssplus)
* A modern introduction to QRSS: [The New Age of QRSS](https://swharden.com/blog/2020-10-03-new-age-of-qrss)
* [FSKview](https://swharden.com/software/FSKview) is a new QRSS and WSPR Spectrogram Viewer for Windows

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.

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 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, τcell): 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, τclamp): 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 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 if you wish to simulate them yourself.

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.

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.

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.

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.

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 too if you want to change some of the parameters and see how it affects the trace.

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. Here we'll analyze the .raw file generated by the whole-cell circuit model above.

# 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

# 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.

⚠️ 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...

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

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:

It = Id * exp(-t / τclamp) + Iss

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:

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 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.

It = m * exp(-t / tau) + b

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.

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.

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)
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%

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 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) we find that cell voltage lags slightly behind the clamp voltage.

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).

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.

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.

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%

Resources

Markdown source code last modified on January 18th, 2021
---
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="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
Pages