Object Type: hsolve Description: An object used for implementing the Hines method for the implicit solution of dendritic trees, as well as a set of routines that maximize speed. This allows faster, more stable numerical integration methods to be used with GENESIS, particularly when there are many compartments in a cell. Author: U. S. Bhalla and E. De Schutter Caltech 91-93, E. De Schutter BBF-UIA 94-99. ------------------------------------------------------------------------------ ELEMENT PARAMETERS DataStructure: hsolve_type [in src/hines/hines_struct.h] Size: 552 bytes Fields: path wildcard list of compartments to be integrated; specifies all compartments belonging to cell which hsolve element will deal with. The default is ./##[][TYPE=compartment] chanmode flag that controls channel computations and chip-array use. 0 = default. It is the most general, uses least memory and is slowest (no chip-array). It takes over the actions of compartments only, computing all other object types as before. As a consequence, all computed fields of the original elements are updated, and all user-setable fields may be set, just as before the element was taken over by the hsolve element. This means that you can add and delete outgoing messages to compartments or other elements whenever you like and easily change parameters during the course of the simulation, making this the easiest mode to use. It is also the most compatible mode of operation, and is guaranteed to work with any future new object type. One significant limitation is that you cannot add or delete AXIAL, RAXIAL or CHANNEL messages, once the hsolve element has been created. In order to calculate the compartment Im field, the compute_Im flag must be set. 1 = like chanmode 0, but optimizes calculations for tabchannel equations. This will be faster than chanmode 0 for models containing tabchannels and tab2Dchannels. 2 = assumes integer exponents (maximum = 6) for tabchannel and tab2Dchannel gate variables. Chanmodes 2-5 are the fastest because of the large 'chip-array'. Under these moses, you can no longer assume that all the fields of the elements that are taken over by hsolve will be updated. Incoming and outgoing messages to and from the disabled elements will work properly, provided that they were added prior to setting up the hsolve element. However, the Vm fields of all compartments are automatically updated under chanmode 2, whether or not there were pre-existing messages to output Vm. Note also that several fields (Gk, Ik, Ek, and Im) are not available for output in chanmodes 2 and 3, so you have to use chanmode 4 or 5 with findsolvefield if you want to output these fields. 3 = as chanmode 2, but Vm fields of compartments are not automatically updated. This mode is faster than chanmode 2, at the expense of somewhat greater memory usage. All disabled elements (including compartments) that had outgoing messages to non-hsolved objects during the SETUP call are updated at a rate determined by outclock. Thus SAVE and PLOT messages will work properly, except with the fields Gk, Ik, Ek, and Im. 4 = as chanmode 3, but more variables (e.g. Gk, Ik, Ek, Im and leak) are stored by hsolve so that they can be accessed by SAVE or PLOT messages (see findsolvefield documentation for more details). The added field, leak, is available for each compartment to give the current flowing through the membrane resistance Rm. Im is calculated directly as the sum of the channel currents and the leakage current. 5 = as chanmode 4, but normalizes the currents and conductances stored in the hsolve givals array (for display purposes only). computeIm flag to determine whether the compartment Im field is calculated in chanmodes 0 and 1. If set, Im is approximated by the sum of axial currents and injected current, as done for non-hsolved compartments and symcompartments. 0 = the default, does not calculate Im. 1 = calculates Im, with a decrease in speed. comptmode flag that selects option for compartment computations 0 = uses least memory, and is slightly slower. However the amount of time spent in this part is usually less than 10% of the total, so this should not matter. 1 = the default, uses a lot of memory and is twice as fast for the actual Hines matrix calculation. This may be only a 5% difference. comptmode 1 is required for chanmodes 2-5. calcmode flag affecting operations for chanmodes 2-5. 0 = no interpolation will be performed in lookup tables. This mode is for backward compatibility with older versions of hsolve (prior to GENESIS 2.1) which did not use interpolation. 1 = default mode with linear interpolation of values in lookup tables. storemode flag, only valid for chanmodes 4 and 5, to allow the output of total currents and conductances. For each type of voltage gated channel in the model, the total currents or conductances are the sum of the corresponding Ik or Gk fields for all compartments where the channel is present. This assumes that these channels have the same name in each compartment. These are stored in an array called itotal. When the hsolve element is set up, a message will be output giving a list of channel names and corresponding itotal indices, if the 'silent' command has previously been given with a negative argument. 0 = no sums are stored. 1 = total currents are stored. 2 = total conductances are stored. no_elminfo = 0/1; flag that controls whether the HPUT and HSET actions will work. Is relevant only for chanmodes 2-5. Default is zero (HPUT and HSET work). Saves memory if non-zero value is used (HPUT, HSET do not work; HRESTORE and HSAVE do). outclock number [0-99] of the clock used for all element updates if chanmode = 3, 4, or 5, will affect SAVE and PLOT messages. The element update routines are not very efficient, so setting them to a slower clock will speed up things. Note however that any object that is not computed by the hines solver and that depends on a message from a hsolve-computed object will be affected. The other fields displayed with the showobject command are NOT to be set by the user. ------------------------------------------------------------------------------ SIMULATION PARAMETERS Function: HinesSolver [in src/hines/hsolve.c] Classes: hsolver Actions: RESET does the standard reset functions, but also for chanmodes 2-5, it will update all parameters in the chip-array (equivalent to a HRESTORE call) and recompute rate factor tables if clocks were changed. DELETE standard actions PROCESS standard actions SETUP setup all internal tables. Must be called before the hsolve element is used. User accessible fields (path, chanmode, etc.) should have been set. Model changes after the SETUP call may have no effect on the computations. DUPLICATE does an efficient duplication of an hsolve element. Use this in a network simulation for identical copies of a neuron. Only the tables containing changing parameters (Vm, etc.) will be duplicated. Syntax: call hsolve1 DUPLICATE hsolve2 duplic_path HPUT updates chip-array. Fields from a single element are put into the chip-array (chanmodes 2 or 3). Do this call after a setfield command on the element. Syntax: call hsolve HPUT element_path HGET updates hsolve-computed element. Computed values are put from the chip-array into a single element (chanmodes 2-5). Do this call before a getfield on the element. Syntax: call hsolve HGET element_path HRESTORE the complete chip-array is updated. Fields from all hsolve-computed elements are put into chip-array (chanmodes 2-5). Do this call after a restore command or after multiple setfield commands. Syntax: call hsolve HRESTORE HSAVE all hsolve-computed elements are updated. Computed values are put from the chip array into all elements (chanmodes 2-5). Do this call before a save command or before multiple getfield commands. Syntax: call hsolve HSAVE Messages: None. Example: // Do all the preparatory grunge work . . // create the cell as an hsolve element readcell test.p /test -hsolve // set the chanmode, comptmode and calcmode (if other than defaults) setfield /test chanmode 2 // set up the arrays and tables for the solver call /test/solve SETUP // use the Crank-Nicholson method for the hsolve element setmethod 11 // It is essential to call reset (or reschedule) after setting // up an hsolver, so that the process list gets updated. reset ------------------------------------------------------------------------------ Notes: The hsolve element completely takes over the calculations for compartment and symcompartment elements and certain other elements specified as part of the cell. Each hsolve element should only solve one cell. The timestep for the integration is determined by the clock assigned to the hsolve element, and the clocks for these elements which are taken over are ignored. When using the Hines solver with a neuron, it is best to think of the entire neuron as a single object since the individual compartments within the cell are no longer responsible for their own computations. This loss of object orientedness is mitigated to some degree by the ability of the solver to transparently interact with elements utilizing other integration schemes. There are 2 basic modes for hsolve operation: without chip-array (chanmode 0 or 1), or with chip-array (chanmodes 2-5). Without chip-array is the most compatible mode, but is the slowest. With the chip-array, hsolve is much faster because the original elements are no longer used, instead all simulation parameters are stored in a huge array (this improves memory access times). Unfortunately you can no longer expect that it will update the fields in computed objects (like Vm in a compartment or Ca in Ca_concen) so that graphic or file output might not work. Vice versa, if you change a field (like inject in compartment) it might not affect the simulation. However, there are methods available to get values in and out of the chip-array, which involve special settings and/or the use of the HPUT, HGET, HSAVE and HRESTORE actions listed above. This is ilustrated in another example below. The documentation for findsolvefield describes another method, introduced in GENESIS 2.1, that also allows access to these values. The example above illustrates the process of setting up the hsolve element. First, one should create the cell as an hsolve element at the root of the cell element tree. This may be done in one step by using the "-hsolve" option with readcell. Note that, starting with GENESIS version 2.2, it is necessary to provide the full path to the cell when using readcell to directly create an hsolve. (e.g. you can't use "readcell test.p test -hsolve", even if "/" is the current working element.) If readcell is not used, then the hsolved cell should be created with statements like create hsolve /test create compartment /test/soma create tabchannel /test/soma/Na_channel ... It is required that only compartments (or symcompartments) be children of the cell, and that channel, concentration, etc. elements should be children or grandchildren of the compartment to which they are attached. Next, set any fields of the hsolve which are needed to specify non-default values of chanmode, comptmode, calcmode, storemode, etc. Then, one needs to call the SETUP action in order to tell the hsolve element to create all the solution arrays and tables. Finally, chose either method 10 (backwards Euler) or 11 (Crank-Nicholson) as the method to be used. All elements that lie outside this tree will continue to be treated by the previous explicit method (typically exponential Euler). The Scripts/examples/hines directory contains an example that which demonstrates both chanmodes 0 and 2. Although it is not presently necessary to create the hsolve element as the root of the element tree, future versions of hsolve may require this. The older method of creating the hsolve is to create the cell as a neutral element, create the hsolve element as a child element, and then set the hsolve path field to indicate which elements will be taken over. Typically, a wildcard path is used to refer to all compartment elements in the simulation. This includes all hsolvable sublelements of these compartments. For example, // readcell will place the compartments below the neutral '/test' readcell test.p /test // create the hines solver element below the cell '/test' create hsolve /test/solve // Specify the path for the solver setfield /test/solve path /test/##[][TYPE=compartment] RESTRICTIONS It is important to be aware of some of the restrictions imposed by the use of hsolve: Only the objects compartment, symcompartment, tabchannel, tab2Dchannel, tabcurrent, spikegen, Ca_concen, concpool, nernst, Mg_block, ghk, taupump, mmpump, hillpump, difshell, fixbuffer, difbuffer, dif2buffer, fura2, synchan, and the oldconn library channels channelC2 and channelC3 are handled by hsolve. If your simulation uses these listed objects, the use of hsolve will increase speed significantly. Note that as of GENESIS version 2.2, symcompartment objects are handled by hsolve. An important restriction introduced in GENESIS version 2.1 is that, for chanmodes 2-5, the element tree of your cell (or other element tree to be taken over by hsolve) must not contain any non-hsolvable elements other than neutral elements. For existing simulations which violate this restriction, it will be easiest to switch to chanmode 0 or 1. It is not guaranteed that messages from unlisted objects to listed objects will work, though solve should alert you about any such incompatibilities. Most incoming messages will work if you use chanmodes 2-5, however you will receive a warning if the source of the message is not hsolved. This is to let you know that the source element is being solved by a less accurate explicit method. Messages from listed to unlisted objects will only work for certain chanmode settings. hsolve does not keep track of changes made to the model after the hsolve SETUP command (particularly the following commands: create, copy, delete, addmsg, deletemsg). In chanmode 0 or 1 this only affects the listed objects and messages between them; in chanmodes 2-5 this affects the listed objects and ALL messages going in or out from them. The hsolver should be used in these modes only for finished models (you can use it for parameter searches); do not use it if you are still constructing and testing a model. hsolve may miscalculate if you change clocks used by tabchannels, tab2Dchannels, synchans, or channelC2/C3 objects without doing a reset (only in chanmodes 2-5). The Hines solver utilizes a considerable amount of memory. Memory use increases as the result of comptmode + chanmode - no_elminfo flags. In particular, going from chanmode 0 or 1 to chanmode 2, 3, or 4 causes a big jump in memory usage. Approximately 100 bytes are required per compartment if the speed-optimized version of the hsolve element is used. Most of this storage is for tables describing the sequence of calculations. These tables can, however, be shared between cells which are identical in topology and differ only in parameter values. FURTHER EXAMPLES * Copying cells with the hsolve element When a cell using the hsolve element is copied, the pointers in the duplicated hsolve element are unchanged. In other words, they refer back to the original cell. One could simply change the path to refer to the current cell and call SETUP again, but that would unnecessarily duplicate a lot of tables. The efficient option would be to use the original tables where they are identical, and construct new ones where they refer to the current cell. This is done using the DUPLICATE sction, which takes the name of the new hsolve element and a wildcard path duplic_path as arguments. duplic_path points to all compartments that should be taken over by the new solver. For example, // Copy the original cell to /test2 copy /test /test2 // create the hines solver element below the cell '/test' create hsolve /test/solve // Specify the path for the solver and the chanmode setfield /test/solve path /test/##[][TYPE=compartment] chanmode 4 call /test/solve SETUP // Duplicate the hsolver call /test/solve DUPLICATE /test2/solve /test2/##[][TYPE=compartment] reset There are several commands, including createmap and cellsheet, that copy cells. The DUPLICATE command should be issued for the hsolvers on each copy. When using chanmodes 3 or higher, findsolvefield must be used with messages in order to access fields of duplicated hsolves, as shown in the documentation for findsolvefield. * Deleting hsolve elements When an hsolve element is set up it removes the relevant channels and compartments from the list of elements to be processed, by setting a bit (0x100) on the flag field of those elements. These elements are re-enabled and the bit set back to 0 when the hsolve element is deleted. WARNING: in the present implementation, the hsolve element does not check whether it has had duplicates made using the DUPLICATE command. This makes it possible to delete the tables for all of the copies by deleting any one of them, resulting in segmentation violations and similar amusements. In other words, do not delete hsolves on duplicated cells. * Accessing data fields in chanmodes 2-5. A script similar to this example might be used for running a long simulation as a background job, using no graphics. The HPUT and HSAVE actions are used in chanmode 2 or 3 to allow the simulation results and the final state of the simulation to be saved to disk. A disk_out element is used to output the Vm of a soma compartment in a compressed binary format. However, it is equally possible to output Vm from every compartment, or from an array of cells. Later, a script using a disk_in element can display the results using messages from the disk_in to an xgraph or xview. The save command is used here to save field values of all elements except nernst elements and the hsolve element. This will allow the simulation to be resumed from its final state by using the restore command. For more details, see Scripts/examples/XODUS/fileview and the documentation for disk_out, disk_in, save, and restore. /* load scripts and global variables */ include defaults include other_stuff pushe /library make_my_stuff pope /* make the model */ readcell my_cell {my_cellpath} -hsolve /* set the clocks */ setclock 0 1e-5 /* integration */ setclock 1 1e-4 /* output */ /* create the output or graphics elements */ /* note that if you want to output Vm only, but from multiple ** compartments, you should use chanmode == 2; if you want to ** output other fields or Vm from only one compartment, ** use chanmode == 3 */ create disk_out /output/disk useclock /output/disk 1 addmsg {my_cellpath}/soma /output/disk Vm addmsg ... /* create any other elements */ .... /* setup the hines solver */ setfield {my_cellpath} \ comptmode 1 \ chanmode 3 \ outclock 1 call {my_cellpath} SETUP setmethod 11 // Crank-Nicholson integration method /* initialize output */ setfield /output/plot_out filename {filename} initialize 1 setfield /output/plot_out filename {filename} append 1 leave_open 1 reset step 0.10 -t /* do a current injection */ setfield {my_cellpath}/soma inject 1.0e-9 /* in Amps */ /* update chip array */ call {my_cellpath} HPUT {my_cellpath}/soma step 2.00 -t /* save the integration values to disk */ call {my_cellpath} HSAVE save {my_cellpath}/##[][TYPE!=nernst][TYPE!=hsolve] {savename} quit --------------------------------------------------------------------------- DETAILED DESCRIPTION The Hines library provides elements and functions for the efficient implicit solution of the systems of differential equations arising in single-neuron models. The sparse matrix arising from the the branched structure of neurons is ordered by the method described by Michael Hines, which permits it to be solved in order N operations using Gaussian elimination without pivoting. The nonlinear equations resulting from the Hodgkin-Huxley description of ion channels are treated as conditionally linear, and also solved in an efficient, second-order manner. Compatibility with other integration schemes is maintained so that mixed integration schemes are feasible. The element sets up a data structure whereby the only change needed for other identical cells is to reassign pointers to the compartment and channel elements. This avoids having to reallocate and rederive the solving scheme. Solutions are done using gaussian forward and backward elimination without pivoting. The numbering and evaluation sequence ensures that no new off-diagonal terms are formed which might mess up the sparse matrix representation. The sparse matrix is represented by a single 'values' array indexed using the ri and cip arrays. The actual solution system is stored as a linear sequence of operations on the 'values' array, in an attempt to speed things up a bit. The implicit method leads to numerically stable solutions for all time steps. Numerical accuracy, however, imposes limits on the integration time step which must be empirically evaluated for each simulation. Typical values are 100 microsecond time steps for purely passive membranes and 20 microsecond time steps for membranes with active channels. This confers a considerable speed up in detailed neuronal simulations where the presence of electrotonically small compartments (such as dendritic spines) leads to systems of equations which are numerically very stiff. In order to obtain similar accuracy using explicit integration methods, time steps in the nanosecond range may be required. In addition to the substantial speed advantage gained from the use of an implicit solution, the coding of the Hines solver element has been done with a view to computational efficiency. This results in a speed up of around 50% compared to the same model with the same time step, using the conventional explicit Euler integration scheme. Code changes needed for the Hines library This section is probably of interest only to individuals who have developed their own code libraries and who are deeply entangled in GENESIS coding issues. For most users the changes are transparent, once the .simrc and schedule.g file have been updated to reflect the addition of a new library. The incorporation of a completely new integration scheme to the simulator can be expected to have numerous repercussions. Fortunately, the changes needed for the Hines library turned out to be relatively painless. No changes were needed for the code of any of the existing elements. There is a slight rearrangement of the order of fields in the various channel objects, and a new globally #define'd type called CHAN_TYPE has been created so that all these channel objects have the following common fields at identical memory offsets: Ik, Gk and Ek (See src/sim/struct_defs.h). All channels whose structures obey this convention can be used in cells utilizing the hsolve element. A very small number of channel objects (such as channelC) do not share these fields and therefore cannot be part of cell models which use the hines solver. A new flag (0x100) has been employed for the Hines solver, which causes the removal of compartment and tabchannel elements from the action list (the list of elements whose action functions are to be called according to the clocks). This flag is similar to the disable option except that it leaves the children of these elements intact. In the interests of speed (which is, after all, what the Hines method is all about) I have 'unrolled' the sparse matrix solution into a single giant function array. This is done at setup time, so the actual solution does not involve any conditionals and minimises array lookups. This is moderately expensive in terms of memory. The size of the function array allocated is S = FA * sizeof (int) ; FA = 10 + 1.5 * M^2 / N where S is the size of the array, FA is the number of functions allocated, M is the number of non-zero coefficients in the solution matrix and N is the number of compartments. M depends strongly on the branching pattern of the cell. FA is an approximation to, and is greater than F, which is the actual number of locations needed in the function table. For example, in a mitral cell model, which has limited branching: N = 286 ; M = 856 ; FA = 3853 ; S = 15412 ; F = 3710 In a granule cell model with numerous dendritic spines: N = 944 ; M = 2830 ; FA = 12736 ; S = 50944 ; F = 12264 A reasonable approximation is that each compartment requires about 54 bytes in the function table, which is not too bad. References: Hines, M. (1984) Efficient computation of branched nerve equations. Int. J. Bio-Med. Comp. 15: 69-76 Mascagni, M.V. (1989) Methods in Neuronal Modeling Ed: Koch and Segev. Chapter 13. 439-484 See also: setmethod, findsolvefield