1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
|
.. _hoc_linmod:
LinearMechanism
---------------
.. hoc:class:: LinearMechanism
Syntax:
``lm = new LinearMechanism(c, g, y, [y0], b)``
``section lm = new LinearMechanism(c, g, y, [y0], b, x)``
``lm = new LinearMechanism(c, g, y, [y0], b, sl, xvec, [layervec])``
``lm = new LinearMechanism(pycallable, c, g, y, ...)``
Description:
Adds linear equations to the tree matrix current balance equations.
I.e. the equations are solved
simultaneously with the current balance equations.
These equations may modify current balance equations and involve
membrane potentials as dependent variables.
The equations added are of the differential-algebraic form
:math:`c \frac{dy}{dt} + g y = b`
with initial conditions specified by the optional y0 vector argument.
c and g must be square matrices of the same rank as the y and b vectors.
The implementation is more efficient if c is a sparse matrix since
at every time step c*y/dt must be computed.
When a LinearMechanism is created, all the potentially non-zero elements
for the c and g matrices must be actually non-zero so that
the mathematical topology of the matrices is known in advance.
After creation, elements can be set to 0 if desired.
The arguments after the b vector specify which voltages and current
balance equations are coupled to this system. The scalar form, x, with
a currently accessed section means that the first equation
is added to the current balance equation at this location and the first
dependent variable is a copy of the membrane potential. If the
system is coupled to more than one location, then sl must be a SectionList
and xvec a Vector of relative positions (0 ... 1) specifying the
locations. In this case, the first xvec.size equations are added to the
corresponding current balance equations and the first xvec.size dependent
y variables are copies of the membrane potentials at this location.
If the optional layervec argument is present then the values must be
0, 1, or 2 (or up to however many layers are defined in :file:`src/nrnoc/options.h`)
0 refers to the internal potential (equal to the membrane potential when
the extracellular mechanism is not inserted), and higher numbers refer
to the \ ``vext[layer-1]`` layer (or ground if the extracellular mechanism is
not inserted).
If some y variables correspond to membrane potential, the corresponding
initial values in the y0 vector are ignored and the initial values come
from the values of v during the normal :hoc:func:`finitialize` call. If you change
the value of v after finitialize, then you should also change the
corresponding y values if the linear system involves derivatives of v.
Note that current balance equations of sections when 0 < x < 1 have dimensions
of milliamp/cm2 and positive terms are outward. Thus
c elements involving voltages in mV
have dimensions of 1000 :math:`\mathrm{\mu{}F/cm^2}` (so a value of .001 corresponds to
1 :math:`\mathrm{\mu{}F/cm^2}`), g elements have dimensions of :math:`\mathrm{S/cm^2}`, and b elements have
dimensions of outward current in :math:`\mathrm{milliamp/cm^2}`. The current balance
equations for the zero area nodes at the beginning and end
of a section (x = 0 and x = 1) have terms with the dimensions of
nanoamps. Thus c elements involving voltages in mV have dimensions
of nF and g elements have dimensions of :math:`\mathrm{\mu{}S}`.
The existence of one or more LinearMechanism switches the gaussian elimination
solver to the general sparse linear equation solver written by
Kenneth S. Kundert and available from
http://www.netlib.org/sparse/index.html
Although, even with no added equations, the solving of m*x=b takes more
than twice as long as the original default solver, there is no restriction
to a tree topology.
Example:
.. code-block::
none
load_file("nrngui.hoc")
create soma
soma { insert hh }
//ideal voltage clamp.
objref c, g, y, b, model
c = new Matrix(2,2,2) //sparse - no elements used
g = new Matrix(2,2)
y = new Vector(2) // y.x[1] is injected current
b = new Vector(2)
g.x[0][1] = -1
g.x[1][0] = 1
b.x[1] = 10 // voltage clamp level
soma model = new LinearMechanism(c, g, y, b, .5)
proc advance() {
printf("t=%g v=%g y.x[1]=%g\n", t, soma.v(.5), y.x[1])
fadvance()
}
run()
.. warning::
Does not work with the CVODE integrator but does work with the
differential-algebraic solver IDA. Note that if the standard
run system is loaded, ``cvode_active(1)`` will automatically
choose the correct variable step integrator.
Does not allow changes to coupling locations.
Is not notified when matrices, vectors, or segments it depends on
disappear.
Description (continued):
If the pycallable argument (A Python Callable object) is present
it is called just before the b Vector is used during a simulation. The
callable can change the elements of b and g (but do not introduce new
elements into g) as a function of time and states. It may be useful for
stability and performance to place the linearized part of b into g.
Consider the following pendulum.py with equations
Example:
.. math::
\frac{d\theta}{dt} = \omega
.. math::
\frac{d\omega}{dt} = -\frac{g}{L} \sin(\theta) \text{ with } \frac{g}{L}=1
.. code-block::
python
from neuron import h
from math import sin
h.load_file('nrngui.hoc')
cmat = h.Matrix(2,2,2).ident()
gmat = h.Matrix(2,2,2)
gmat.setval(0,1, -1)
y = h.Vector(2)
y0 = h.Vector(2)
b = h.Vector(2)
def callback():
b.x[1] = -sin(y.x[0])
nlm = h.LinearMechanism(callback, cmat, gmat, y, y0, b)
dummy = h.Section()
trajec = h.Vector()
tvec = h.Vector()
trajec.record(y._ref_x[0])
tvec.record(h._ref_t)
graph = h.Graph()
h.tstop=50
def prun(theta0, omega0):
graph.erase()
y0.x[0] = theta0
y0.x[1] = omega0
h.run()
trajec.line(graph, tvec)
h.dt /= 10
h.cvode.atol(1e-5)
h.cvode_active(1)
prun(0, 1.9999) # 2.0001 will keep it rotating
graph.exec_menu("View = plot")
|
