The voltage jump method

A Mathematica notebook

This Mathematica notebook is an implementation of Equations 1-10 in Häusser M, Roth A (1997) Estimating the time course of the excitatory synaptic conductance in neocortical pyramidal cells using a novel voltage jump method, J Neurosci 17:7606-7625. The input cells should be evaluated in order, starting from a fresh Mathematica session. Some cells may take several minutes to evaluate, depending on the computer system used. This notebook was written for Mathematica version 3.0 or greater.

Authors:	Michael Häusser	m.hausser@ucl.ac.uk
		Arnd Roth		roth@sunny.mpimf-heidelberg.mpg.de
Version 1.1	22.5.1998

<< "Calculus`DiracDelta`"Off[General::"spell"]; Off[General::"spell1"];

Analytical functions describing the charge recovery

Delta function conductance (Eq . 1 - 5)

Vsyn := 
      alpha UnitStep[t - s] 
        (Vcom - Vcom (E^(-((t - s)/tauv))))

g := gHat DiracDelta[t]

Isyn := Vsyn g

Qsyn = 
      (∫_(-∞)^∞) Isyn \[DifferentialD]t

Qsoma = alpha Qsyn

Plot[
      Evaluate[Qsoma /. {alpha -> 0.85, tauv -> 2.0, 
            gHat -> 1.0, Vcom -> -1.0}], {s, -10, 10}, 
      AxesOrigin -> {-10, 0}, 
      AxesLabel -> {"s", "Qsoma"}, PlotPoints -> 200, 
      PlotStyle -> {{Thickness[0.01]}}]

Single - exponential decay of conductance (Eq . 6 - 7)

g := gBar UnitStep[t] (E^(-(t/taudec)))

myIntegrate[
        f_ UnitStep[t_] UnitStep[t_ - s_], {t_, -∞, 
          ∞}] := 
      If[s ≤ 0, 
        Evaluate[
          Simplify[
            Integrate[f, {t, 0, ∞}, 
              GenerateConditions -> False]]], 
        Evaluate[
          Simplify[
            Integrate[f, {t, s, ∞}, 
              GenerateConditions -> False, Assumptions -> s > 0]]]]

Qsoma = 
      myIntegrate[alpha Isyn, {t, -∞, ∞}]

Plot[
      Evaluate[Qsoma /. {alpha -> 0.85, tauv -> 2.0, 
            gBar -> 1.0, Vcom -> -1.0, taudec -> 1.5}], {
        s, -10, 10}, AxesOrigin -> {-10, 0}, 
      AxesLabel -> {"s", "Qsoma"}, PlotPoints -> 200, 
      PlotStyle -> {{Thickness[0.01]}}]

More general case (Eq . 8 - 10)

Vsyn := 
      alpha Vcom UnitStep[t - s] 
        (av1 + av2 - av1 (E^(-((t - s)/tauv1))) - 
            av2 (E^(-((t - s)/tauv2))))

g := 
      UnitStep[t] 
        (-(gBar1 + gBar2) (E^(-(t/taurise))) + 
            gBar1 (E^(-(t/taudec1))) + gBar2 (E^(-(t/taudec2))))

Qsoma = 
      myIntegrate[alpha Isyn, {t, -∞, ∞}]

Plot[
      Evaluate[Qsoma /. {alpha -> 0.85, av1 -> 0.7, 
            tauv1 -> 2.0, av2 -> 0.3, tauv2 -> 5.0, 
            taurise -> 0.2, gBar1 -> 0.9, taudec1 -> 1.5, 
            gBar2 -> 0.1, taudec2 -> 5.0, Vcom -> -1.0}], {
        s, -10, 10}, AxesOrigin -> {-10, 0}, 
      AxesLabel -> {"s", "Qsoma"}, PlotPoints -> 200, 
      PlotStyle -> {{Thickness[0.01]}}]

Fitting charge recovery data

Simplifying Eq. 10: single-exponential rise and decay of conductance

ClearAttributes[If, HoldRest]Qsoma = 
      Simplify[Qsoma /. {alpha -> 1, av1 -> 1, tauv1 -> tauv, 
            av2 -> 0, gBar1 -> gStar, taudec1 -> taudec, 
            gBar2 -> 0, Vcom -> -1}]SetAttributes[If, HoldRest]

Implementation of fit procedure

Off[FindMinimum::"precw"]; 
    fitFunction = 
      Compile[{tauv, gStar, taurise, taudec, s}, 
        If[s ≤ 0, 
          -((gStar (taudec - taurise) 
                  (taudec (taurise + tauv) + 
                      tauv (taurise + tauv - (E^(s/tauv)) tauv))
                    )/((taudec + tauv) (taurise + tauv))), 
          ((gStar (
                (E^(-(s/taurise))) taudec (taurise^2) + 
                  (E^(-(s/taurise))) (taurise^2) tauv - 
                  (E^(-(s/taudec))) (taudec^2) (taurise + tauv))
                )/((taudec + tauv) (taurise + tauv)))]]; 
    chiSquare[{tv_, gb_, trise_, tdec_, xoffset_, yoffset_}, data_List] := 
      Plus@@(
            (((Last[#1] - 
                    fitFunction[tv, gb, trise, tdec, First[#1] - xoffset] - 
                    yoffset))^2)&)/@datatrisewindow[trise_] := 
      Which[trise < 0.1, 1000.  (((trise - 0.1))^2), trise > 0.6, 
        1000.  (((trise - 0.6))^2), True, 0. ]chargeRecoveryFit[
        data_List, {gtv_, gbg_, gtrise_, gtdec_, gxoffset_, gyoffset_}] := 
      FindMinimum[
        chiSquare[{tv, gb, trise, tdec, xoffset, yoffset}, data] + 
          trisewindow[trise], {tv, {gtv, 1.05 gtv}, 0. , 1000. }, {
          gb, {gbg, 1.05 gbg}, 0. , 1000. }, {tdec, {gtdec, 1.05 gtdec}, 
          0. , 1000. }, {xoffset, {gxoffset, 0.05 + gxoffset}, -100. , 
          100. }, {yoffset, {gyoffset, 0.05 + gyoffset}, -100. , 100. }, {
          trise, {gtrise, 1.05 gtrise}, 0.02, 1000. }, 
        PrecisionGoal -> 10, AccuracyGoal -> 10, 
        MaxIterations -> 10000]

Example fit

data = 
      Table[{s, 
          Qsoma /. {tauv -> 2.0, gStar -> (0.85^2), 
              taurise -> 0.5, taudec -> 1.5}}, {s, -10, 10, 
          0.5}]

dataPlot = 
      ListPlot[data, AxesOrigin -> {-10, 0}, 
        AxesLabel -> {"s", "Qsoma"}, 
        PlotStyle -> {PointSize[0.03]}]

fit = chargeRecoveryFit[data, {2.2, 0.8, 0.35, 1.4, 0.1, 0.01}]

Off[CompiledFunction::"cfr"]fitPlot = 
      Plot[Evaluate[
          fitFunction[tv, gb, trise, tdec, s - xoffset] + yoffset /. 
            Last[fit]], {s, -10, 10}, Compiled -> False, 
        AxesOrigin -> {-10, 0}, 
        AxesLabel -> {"s", "Qsoma"}, PlotPoints -> 200, 
        PlotStyle -> {{Thickness[0.01]}}]On[CompiledFunction::"cfr"]

Show[dataPlot, fitPlot]