Page 1 of 1

error in saving a voltage trace into a figure

Posted: Tue Sep 24, 2019 11:20 am
by alexandrapierri
Hello developers

I am trying to do the simple task, save a figure using one of the python templates found on the NEURON Documentation: initonerun.py.
The script plots the Graph just fine using one of my own cell templates (example.hoc), but it doesn't save it correctly. The problem is in the last six lines of the code from "#issue starts here" to the end, which are my own addition.

I wrote these lines to save my Graph into a figure (the total voltage trace as a function of time) consulting the Programmer's Reference and the forum, but in doing so I generate an empty figure with wrong ranges for both the voltage and the time axis. What am I doing wrong ?

thanks for all help,
Alexandra

Code: Select all

"""
Executes one simulation with specified stimulus.
Displays response and reports spike frequency
"""
from neuron import h, gui
import matplotlib.pyplot as plt
import numpy as np

# simulation parameter
h.tstop = 500 # ms, more than long enough for 15 spikes at ISI = 25 ms

# model specification
#from cell import Cell

h.load_file('example.hoc') 

cell = h.example() 
soma_loc=cell.soma[0]
dend_loc=cell.dend[1]
# instrumentation

# experimental manipulations
pre=h.Section()
stim = h.IClamp(soma_loc(0.5,sec=pre))
stim.delay = 10 # ms
stim.dur = 20e9
stim.amp = 0.5 # nA

# data recording and analysis
# count only those spikes that get to distal end of dend
nc = h.NetCon(dend_loc(1)._ref_v, None, sec=dend_loc)
nc.threshold = -10 # mV
spvec = h.Vector()
nc.record(spvec)

NSETTLE = 5 # ignore the first NSETTLE ISI (allow freq to stablize)
NINVL = 10 # num ISI from which frequency will be calculated
NMIN = NSETTLE + NINVL # ignore recordings with fewer than this num of ISIs

def get_frequency():
    nspikes = spvec.size()
    if nspikes > NMIN:
        t2 = spvec[-1] # last spike
        t1 = spvec[-(1 + NINVL)] # NINVL prior to last spike
        return NINVL * 1.e3 / (t2 - t1)
    else:
        return 0

# Simulation control and reporting of results

def onerun(amp):
    # amp = amplitude of stimulus current
    g = h.Graph(0)
    g.size(0, 500, -80, 40)
    g.view(0, -80, 500, 120, 2, 105, 300.48, 200.32)
    # update graph throughout the simulation
    h.graphList[0].append(g)
    # plot v at distal end of dend
    g.addvar('dend(1).v', dend_loc(1)._ref_v)
    stim.amp = amp
    h.run()
    freq = get_frequency()
    print('stimulus %g frequency %g' % (amp, freq))
    g.save_name("graphList.")
onerun(0.15)

vvecA = h.Vector().record(dend_loc(1)._ref_v)
tvec = h.Vector().record(h._ref_t)

#Issue starts here 
fig = plt.figure()
plt.plot(tvec, vvecA, label="it_worked!")
plt.xlabel('t (ms)')
plt.ylabel('V$_m$ (mV)')
plt.legend(frameon=False)
fig.savefig("trial.png")

Re: error in saving a voltage trace into a figure

Posted: Tue Sep 24, 2019 12:14 pm
by vogdb
Hi! Does it work fine if you plot directly instead of saving figure? What would happen if you replace

Code: Select all

fig.savefig("trial.png")
with

Code: Select all

plt.show()
?

Re: error in saving a voltage trace into a figure

Posted: Tue Sep 24, 2019 12:19 pm
by ted
Good question. Many people make similar mistakes. Here's a chance to learn a lot.
What am I doing wrong?
Misdiagnosing the cause of the problem because of failure to localize the problem.

In medicine, the first step in diagnosis is to localize the cause of the problem. This is also true for diagnosing sofware bugs. The symptom you report is an empty graph. Are you quite sure that the "issue" starts where you think it does? Maybe the code below
#Issue starts here
is correct--it simply reveals the symptom--and the problem originates in earlier lines. Have you confirmed that vvecA and tvec even contain data? If they don't, why don't they? Hint: Read the Programmer's Reference documentation of the Vector class's record() method, then think about

Code: Select all

def onerun(amp):
  . . .
  h.run()
  . . .

onerun(0.15)

vvecA = h.Vector().record(dend_loc(1)._ref_v)
tvec = h.Vector().record(h._ref_t)
After you fix that, read the documentation of the IClamp class, then explain why it does or doesn't make sense to write
pre=h.Section()
stim = h.IClamp(soma_loc(0.5,sec=pre))
and revise these statements as appropriate. This isn't causing any symptoms, but it's not correct, and in a future version of NEURON it might cause problems.