Skip to content

Commit

Permalink
wip(gfluct): conversion complete, testing in progress
Browse files Browse the repository at this point in the history
  • Loading branch information
@sanjayankur31
sanjayankur31 committed Aug 29, 2023
1 parent aaad187 commit dcc70f7
Show file tree
Hide file tree
Showing 2 changed files with 150 additions and 71 deletions.
151 changes: 136 additions & 15 deletions NeuroML2/Gfluct.nml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,87 @@

<!-- Note that order of different LEMS types matters for validation, since the schema defines what order they should appear in -->

<!-- Reference:
<!--
Fluctuating conductance model for synaptic bombardment
======================================================
THEORY
Synaptic bombardment is represented by a stochastic model containing
two fluctuating conductances g_e(t) and g_i(t) descibed by:
Isyn = g_e(t) * [V - E_e] + g_i(t) * [V - E_i]
d g_e / dt = -(g_e - g_e0) / tau_e + sqrt(D_e) * Ft
d g_i / dt = -(g_i - g_i0) / tau_i + sqrt(D_i) * Ft
where E_e, E_i are the reversal potentials, g_e0, g_i0 are the average
conductances, tau_e, tau_i are time constants, D_e, D_i are noise diffusion
coefficients and Ft is a gaussian white noise of unit standard deviation.
g_e and g_i are described by an Ornstein-Uhlenbeck (OU) stochastic process
where tau_e and tau_i represent the "correlation" (if tau_e and tau_i are
zero, g_e and g_i are white noise). The estimation of OU parameters can
be made from the power spectrum:
S(w) = 2 * D * tau^2 / (1 + w^2 * tau^2)
and the diffusion coeffient D is estimated from the variance:
D = 2 * sigma^2 / tau
NUMERICAL RESOLUTION
The numerical scheme for integration of OU processes takes advantage
of the fact that these processes are gaussian, which led to an exact
update rule independent of the time step dt (see Gillespie DT, Am J Phys
64: 225, 1996):
x(t+dt) = x(t) * exp(-dt/tau) + A * N(0,1)
where A = sqrt( D*tau/2 * (1-exp(-2*dt/tau)) ) and N(0,1) is a normal
random number (avg=0, sigma=1)
IMPLEMENTATION
This mechanism is implemented as a nonspecific current defined as a
point process.
PARAMETERS
The mechanism takes the following parameters:
E_e = 0 (mV) : reversal potential of excitatory conductance
E_i = -75 (mV) : reversal potential of inhibitory conductance
g_e0 = 0.0121 (umho) : average excitatory conductance
g_i0 = 0.0573 (umho) : average inhibitory conductance
std_e = 0.0030 (umho) : standard dev of excitatory conductance
std_i = 0.0066 (umho) : standard dev of inhibitory conductance
tau_e = 2.728 (ms) : time constant of excitatory conductance
tau_i = 10.49 (ms) : time constant of inhibitory conductance
Gfluct2: conductance cannot be negative
REFERENCE
Destexhe, A., Rudolph, M., Fellous, J-M. and Sejnowski, T.J.
Fluctuating synaptic conductances recreate in-vivo-like activity in
neocortical neurons. Neuroscience 107: 13-24 (2001).
(electronic copy available at http://cns.iaf.cnrs-gif.fr)
A. Destexhe, 1999
Reference:
https://doi.org/10.1016/S0306-4522(01)00344-X
Destexhe, A.; Rudolph, M.; Fellous, J.-M. & Sejnowski, T. J. Fluctuating synaptic conductances recreate in vivo-like activity in neocortical neurons Neuroscience, 2001, 107, 13-24
Original mod source: https://modeldb.science/8115
-->
Expand All @@ -31,20 +109,39 @@
<Parameter name="tau_e" dimension="time" description="Time constant of excitatory conductance"/>
<Parameter name="tau_i" dimension="time" description="Time constant of inhibitory conductance"/>

<Constant name="zeroms" value="0 ms" dimension="time" />

<EventPort name="in" direction="in" description="Note this is not used here. Will be removed in future"/>

<DerivedParameter name="exp_e" dimension="none" value="exp(-dt/tau_e)" />
<DerivedParameter name="amp_e" dimension="conductance" value="std_e * sqrt((1 - exp(-2 * (dt/tau_e))))" />
<DerivedParameter name="exp_i" dimension="none" value="exp(-dt/tau_i)" />
<DerivedParameter name="amp_i" dimension="conductance" value="std_i * sqrt((1 - exp(-2 * (dt/tau_i))))" />

<Exposure name="g_e" dimension="conductance" />
<Exposure name="g_i" dimension="conductance" />
<Exposure name="g_e1" dimension="conductance" />
<Exposure name="g_i1" dimension="conductance" />
<Exposure name="i" dimension="current" />

<Dynamics>
<!-- total conductances -->
<StateVariable name="g_e" exposure="g_e" dimension="conductance" />
<StateVariable name="g_i" exposure="g_i" dimension="conductance" />
<!-- fluctuating conductances -->
<StateVariable name="g_e1" dimension="conductance" />
<StateVariable name="g_i1" dimension="conductance" />
<StateVariable name="g_e1" exposure="g_e1" dimension="conductance" />
<StateVariable name="g_i1" exposure="g_i1" dimension="conductance" />

<!-- track old value of g_e1: LEMS doesn't seem to like it if the same variable is used -->
<StateVariable name="g_e1_old" dimension="conductance" />
<StateVariable name="g_i1_old" dimension="conductance" />

<StateVariable name="i" exposure="i" dimension="current" />

<StateVariable name="r1" dimension="none" />
<StateVariable name="r2" dimension="none" />
<StateVariable name="randn" dimension="none" />

<OnEvent port="in"><!--TODO: remove, see above...
<StateAssignment variable="i" value="0"/>-->
</OnEvent>
Expand All @@ -55,37 +152,59 @@
<StateAssignment variable="g_e" value="0"/>
<StateAssignment variable="g_i" value="0"/>
<StateAssignment variable="g_e1" value="0" />
<StateAssignment variable="g_i1" value="0" />
<StateAssignment variable="g_e1_old" value="0" />
<StateAssignment variable="g_i1_old" value="0" />
</OnStart>

<!-- before start -->
<OnCondition test="t .lt. start">
<StateAssignment variable="i" value="0"/>
<StateAssignment variable="g_e" value="0"/>
<StateAssignment variable="g_i" value="0"/>
<StateAssignment variable="g_e1" value="0" />
<StateAssignment variable="g_i1" value="0" />
<StateAssignment variable="g_e1_old" value="0" />
<StateAssignment variable="g_i1_old" value="0" />
</OnCondition>

<!-- generate gaussian random number -->
<OnCondition test="t .geq. start .and. t.lt. stop">
<StateAssignment variable="r1" value="random(1)"/>
<StateAssignment variable="r2" value="random(1)"/>
<StateAssignment variable="randn" value="sqrt(-2*log(r1))*cos(2*3.14159265359*r2)"/>
</OnCondition>

<!-- oup()-->
<OnCondition test="tau_e .neq. 0" >
<!-- let exp_e and std_e be computed for each step: we can't have conditionals inside OnStart, although I could use H etc. to hack it -->
<StateAssignment variable="g_e" value="exp(-dt/tau_e) * g_e1 + ((std_e * sqrt(1 - exp(-2 * dt/tau_e))) * random(1))" />
<OnCondition test="(t .geq. start .and. t .lt. stop) .and. (tau_e .neq. zeroms)" >
<StateAssignment variable="g_e1" value="(exp_e * g_e1_old) + (amp_e * randn)" />
<StateAssignment variable="g_e1_old" value="g_e1" />
</OnCondition>

<OnCondition test="tau_i .neq. 0" >
<StateAssignment variable="g_i" value="exp(-dt/tau_i) * g_i1 + ((std_i * sqrt(1 - exp(-2 * dt/tau_i))) * random(1))" />
<OnCondition test="(t .geq. start .and. t .lt. stop) .and. (tau_i .neq. zeroms)" >
<StateAssignment variable="g_i1" value="(exp_i * g_i1_old) + (amp_i * randn)" />
<StateAssignment variable="g_i1_old" value="g_i1" />
</OnCondition>

<!-- original mod says g_e = std_e * grand() but that would immediately get overwritten in the next bit, so we think it's the fluctuating current that is calculated here, g_e1 -->
<OnCondition test="tau_e .eq. 0" >
<StateAssignment variable="g_e1" value="std_e * random(1)" />
<OnCondition test="(t .geq. start .and. t .lt. stop) .and. (tau_e .eq. zeroms)" >
<StateAssignment variable="g_e1" value="std_e * randn" />
</OnCondition>
<OnCondition test="tau_i .eq. 0" >
<StateAssignment variable="g_i1" value="std_i * random(1)" />
<OnCondition test="(t .geq. start .and. t .lt. stop) .and. (tau_i .eq. zeroms)" >
<StateAssignment variable="g_i1" value="std_i * randn" />
</OnCondition>

<OnCondition test="t .geq. start .and. t .lt. stop">
<StateAssignment variable="g_e" value="g_e0 + g_e1" />
<StateAssignment variable="g_i" value="g_i0 + g_i1" />
</OnCondition>

<OnCondition test="g_e .lt. 0">
<StateAssignment variable="g_e" value="0" />
</OnCondition>
<OnCondition test="g_i .lt. 0">
<StateAssignment variable="g_i" value="0" />
</OnCondition>

<OnCondition test="t .geq. start .and. t .lt. stop">
<StateAssignment variable="i" value="-1 * g_e * (v - E_e) + g_i * (v - E_i)" />
</OnCondition>

Expand All @@ -94,6 +213,8 @@
<StateAssignment variable="i" value="0"/>
<StateAssignment variable="g_e" value="0"/>
<StateAssignment variable="g_i" value="0"/>
<StateAssignment variable="g_e1" value="0" />
<StateAssignment variable="g_i1" value="0" />
</OnCondition>
</Dynamics>

Expand Down
70 changes: 14 additions & 56 deletions NeuroML2/LEMS_Gfluct_test.xml
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,5 @@
<Lems>


<!-- Specify which component to run -->
<Target component="simnet1"/>

Expand All @@ -18,76 +17,35 @@
tau_syn_E="0.5" tau_syn_I="0.5" v_init="-65" v_reset="-65.0" v_rest="-65.0" v_thresh="-55.0"/>


<pulseGenerator id="pulseGen0" delay="0ms" duration="30000ms" amplitude="0.15 nA" />

<Gfluct id="noisyCurrentSource1" start="0ms" stop="30000ms" dt="0.05ms" E_e="0mV" E_i="-75mV" g_e0="0.0121 uS" g_i0="0.573 uS" std_e="0.003 uS" std_i="0.0066 uS" tau_e="2.728 ms" tau_i="10.49 ms"/>
<Gfluct id="noisyCurrentSource2" start="0ms" stop="30000ms" dt="0.05ms" E_e="0mV" E_i="-75mV" g_e0="0.0121 uS" g_i0="0.573 uS" std_e="0.003 uS" std_i="0.0066 uS" tau_e="2.728 ms" tau_i="10.49 ms"/>
<Gfluct id="noisyCurrentSource3" start="0ms" stop="30000ms" dt="0.05ms" E_e="0mV" E_i="-75mV" g_e0="0.0121 uS" g_i0="0.573 uS" std_e="0.003 uS" std_i="0.0066 uS" tau_e="2.728 ms" tau_i="10.49 ms"/>
<Gfluct id="noisyCurrentSource1" start="1000ms" stop="2000ms" dt="0.5ms" E_e="0mV" E_i="-75mV" g_e0="0.0121 uS" g_i0="0.0573 uS" std_e="0.012 uS" std_i="0.0264 uS" tau_e="2.728 ms" tau_i="10.49 ms"/>

<network id="net1">
<population component="IF_curr_exp" id="Pop0" type="populationList" size="3">
<population component="IF_curr_exp" id="Pop0" type="populationList" size="1">
<instance id="0">
<location x="120" y="230" z="567"/>
</instance>
<instance id="1">
<location x="270" y="450" z="56"/>
</instance>
<instance id="2">
<location x="54" y="234" z="89"/>
</instance>
<instance id="3">
<location x="1" y="2" z="3"/>
</instance>
</population>


<inputList id="stimInput" component="pulseGen0" population="Pop0">
<inputList id="Gfluct0" component="noisyCurrentSource1" population="Pop0">
<input id="0" target="../Pop0/0/IF_curr_exp" destination="synapses"/>
</inputList>

<inputList id="noisy1" component="noisyCurrentSource1" population="Pop0">
<input id="0" target="../Pop0/1/IF_curr_exp" destination="synapses"/>
</inputList>

<inputList id="noisy2" component="noisyCurrentSource2" population="Pop0">
<input id="0" target="../Pop0/2/IF_curr_exp" destination="synapses"/>
</inputList>

<inputList id="noisy3" component="noisyCurrentSource3" population="Pop0">
<input id="0" target="../Pop0/3/IF_curr_exp" destination="synapses"/>
</inputList>

</network>

<Simulation id="simnet1" length="30000ms" step="0.05ms" target="net1">

<Display id="display_voltages" title="Voltages" timeScale="1ms" xmin="-2.0" xmax="10020.0" ymin="-68" ymax="-47">
<Line id="Pop0/0: Vm" quantity="Pop0/0/IF_curr_exp/v" scale="1mV" color="#0000ff" timeScale="1ms"/>
<Line id="Pop0/1: Vm" quantity="Pop0/1/IF_curr_exp/v" scale="1mV" color="#00ff00" timeScale="1ms"/>
<Line id="Pop0/2: Vm" quantity="Pop0/2/IF_curr_exp/v" scale="1mV" color="#ff0000" timeScale="1ms"/>
<Line id="Pop0/3: Vm" quantity="Pop0/3/IF_curr_exp/v" scale="1mV" color="#00ffff" timeScale="1ms"/>
</Display>

<Display id="display_currents" title="Voltages" timeScale="1ms" xmin="-2.0" xmax="10020.0" ymin="-0.05" ymax="0.6">
<Line id="Pop0/0: i" quantity="Pop0/0/IF_curr_exp/pulseGen0/i" scale="1nA" color="#0000ff" timeScale="1ms"/>
<Line id="Pop0/1: i" quantity="Pop0/1/IF_curr_exp/noisyCurrentSource1/i" scale="1nA" color="#00ff00" timeScale="1ms"/>
<Line id="Pop0/2: i" quantity="Pop0/2/IF_curr_exp/noisyCurrentSource2/i" scale="1nA" color="#ff0000" timeScale="1ms"/>
<Line id="Pop0/3: i" quantity="Pop0/3/IF_curr_exp/noisyCurrentSource3/i" scale="1nA" color="#00ffff" timeScale="1ms"/>
</Display>

<OutputFile id="Volts_file" fileName="v_ou.dat">
<OutputColumn id="v0" quantity="Pop0/0/IF_curr_exp/v"/>
<OutputColumn id="v1" quantity="Pop0/1/IF_curr_exp/v"/>
<OutputColumn id="v2" quantity="Pop0/2/IF_curr_exp/v"/>
<OutputColumn id="v3" quantity="Pop0/3/IF_curr_exp/v"/>
</OutputFile>
<Simulation id="simnet1" length="3000ms" step="0.5ms" target="net1">

<OutputFile id="Currents_file" fileName="i_ou.dat">
<OutputColumn id="i0" quantity="Pop0/0/IF_curr_exp/pulseGen0/i"/>
<OutputColumn id="i1" quantity="Pop0/1/IF_curr_exp/noisyCurrentSource1/i"/>
<OutputColumn id="i2" quantity="Pop0/2/IF_curr_exp/noisyCurrentSource2/i"/>
<OutputColumn id="i3" quantity="Pop0/3/IF_curr_exp/noisyCurrentSource3/i"/>
<!-- works for jLEMS -->
<OutputFile id="Currents_file" fileName="i_gfluct.dat">
<OutputColumn id="i" quantity="Pop0/0/IF_curr_exp/synapses:noisyCurrentSource1:0/i"/>
</OutputFile>
<OutputFile id="Conductance_file" fileName="g_gfluct.dat">
<OutputColumn id="g_e" quantity="Pop0/0/IF_curr_exp/synapses:noisyCurrentSource1:0/g_e"/>
<OutputColumn id="g_e1" quantity="Pop0/0/IF_curr_exp/synapses:noisyCurrentSource1:0/g_e1"/>
<OutputColumn id="g_i" quantity="Pop0/0/IF_curr_exp/synapses:noisyCurrentSource1:0/g_i"/>
<OutputColumn id="g_i1" quantity="Pop0/0/IF_curr_exp/synapses:noisyCurrentSource1:0/g_i1"/>
</OutputFile>
<!-- bug: input properties cannot currently be recorded using neuron -->

</Simulation>

Expand Down

0 comments on commit dcc70f7

Please sign in to comment.