Initialization/nseg problem

jaambros · Post by **jaambros** » Fri Jul 14, 2006 11:41 pm

I have a custom init routine that uses the d-lambda rule to set nseg (taking care to always use odd values) and the example in the Neuron book using large dt to find the steady state (see below)

The problem is that repeated executions do not yield the same results.
Eg, in the call sequence init(0.1), init(0.5), init(0.1), the first init(0.1) diverges while the second converges.

It appears that nseg changes can be quite disruptive.
Any suggestions on where I went wrong?

Thanks in advance,
Jose Ambros-Ingerson

Code: Select all

proc init() {   
  mulfit_cell_init( $1 )
  init_steady_state()
}

proc mulfit_cell_init() { 
  forall { // set for lambda_f call
    Ra		= G_Ra
    cm 		= G_cm
  }
  nseg_tot = 0
  forall { 
    nseg = int((L/($1 *lambda_f(100))+.9)/2)*2 + 1 
    nseg_tot += nseg
  }
  printf( "lambda-d %g nseg_tot %d\n", $1, nseg_tot )

  forall {
    e_pas 	= G_e_pas
    g_pas	= 1 / ( G_Rm * 1000 )
    
    gbar_Naf_i0	= Gbar_Naf
    gbar_Nap_i0	= Gbar_Nap
    gbar_KDR_i0 = Gbar_KDR
    gbar_h_i0 = Gbar_h_soma 
    gbar_KA_i0 = Gbar_KA_soma
  }
}

proc init_steady_state(){ local dtsav, was_active
  dtsav	= dt				// record dt and cvode state
  was_active  = cvode.active()
  cvode.active( 0 )
  finitialize( -82.34 )
  fcurrent()
  t     = -3e9
  dt	= 1e8
  
  while( t < 0 ) { 
    if( iss_idebug ) printf( "s1 %g %g\n", t, v )
    fadvance()
  }
  
  cvode.active( was_active )		// restore cvode and dt state
  dt  = dtsav
  t = 0
  if( cvode.active() ) {
    cvode.re_init()
  } else {
    fcurrent()
  }
  frecord_init()
}
... load cell and set parameter values (G_Ra, G_cm, etc)...
init( 0.1 )
init( 0.5 )
init( 0.1 )

Output was

Code: Select all

NEURON -- Version 5.8 2005-10-7 13:46:29 Main (85)
by John W. Moore, Michael Hines, and Ted Carnevale
Duke and Yale University -- Copyright 1984-2005

	1 
	1 
nrn_load_dll /axon/d2/cNeuro/lib/nrn/NMOD/ca1n1-mod/powerpc/.libs/libnrnmech.so
	80 
Additional mechanisms from files
 KA_i0.mod KA_n1.mod KDR_b0.mod KDR_b1.mod KDR_i0.mod KDR_mig.mod Ktf_p1.mod Naf_b1.mod Naf_i0.mod Nap_i0.mod fK_DR_n1.mod fNa_n1.mod h_i0.mod h_n1.mod morpho.mod pNa_n1.mod
loading membrane mechanisms from /axon/d2/cNeuro/lib/nrn/NMOD/ca1n1-mod/powerpc/.libs/libnrnmech.so
	1 
	1 
	1 
	1 
	0.001 
nfeval 34681
	13 
nrniv: unable to open font "*Arial*bold*--12*", using "fixed"
	1 
lambda-d 0.1 nseg_tot 717
s1 -3e+09 -82.34
s1 -2.9e+09 -76.194
s1 -2.8e+09 -73.2505
s1 -2.7e+09 -70.2672
s1 -2.6e+09 -65.9441
s1 -2.5e+09 -56.5226
s1 -2.4e+09 -22.7331
s1 -2.3e+09 -71.4005
s1 -2.2e+09 -60.302
s1 -2.1e+09 -38.9746
s1 -2e+09 -35.0278
s1 -1.9e+09 -50.1262
s1 -1.8e+09 -17.6355
s1 -1.7e+09 -74.357
s1 -1.6e+09 -43.5605
s1 -1.5e+09 -17.8408
s1 -1.4e+09 -74.9825
s1 -1.3e+09 -40.4702
s1 -1.2e+09 -30.904
s1 -1.1e+09 -56.6353
s1 -1e+09 -33.6248
s1 -9e+08 -49.1995
s1 -8e+08 -19.2721
s1 -7e+08 -70.9034
s1 -6e+08 -42.8433
s1 -5e+08 -19.5439
s1 -4e+08 -73.8706
s1 -3e+08 -51.5855
s1 -2e+08 -9.08034
s1 -1e+08 -82.9557
lambda-d 0.5 nseg_tot 251
s1 -3e+09 -82.34
s1 -2.9e+09 -81.2169
s1 -2.8e+09 -81.0766
s1 -2.7e+09 -81.0331
s1 -2.6e+09 -81.0187
s1 -2.5e+09 -81.0134
s1 -2.4e+09 -81.0115
s1 -2.3e+09 -81.0107
s1 -2.2e+09 -81.0104
s1 -2.1e+09 -81.0103
s1 -2e+09 -81.0102
s1 -1.9e+09 -81.0102
s1 -1.8e+09 -81.0102
s1 -1.7e+09 -81.0102
s1 -1.6e+09 -81.0102
s1 -1.5e+09 -81.0102
s1 -1.4e+09 -81.0102
s1 -1.3e+09 -81.0102
s1 -1.2e+09 -81.0102
s1 -1.1e+09 -81.0102
s1 -1e+09 -81.0102
s1 -9e+08 -81.0102
s1 -8e+08 -81.0102
s1 -7e+08 -81.0102
s1 -6e+08 -81.0102
s1 -5e+08 -81.0102
s1 -4e+08 -81.0102
s1 -3e+08 -81.0102
s1 -2e+08 -81.0102
s1 -1e+08 -81.0102
lambda-d 0.1 nseg_tot 717
s1 -3e+09 -82.34
s1 -2.9e+09 -81.1348
s1 -2.8e+09 -80.9879
s1 -2.7e+09 -80.9433
s1 -2.6e+09 -80.9291
s1 -2.5e+09 -80.9242
s1 -2.4e+09 -80.9225
s1 -2.3e+09 -80.9218
s1 -2.2e+09 -80.9216
s1 -2.1e+09 -80.9215
s1 -2e+09 -80.9215
s1 -1.9e+09 -80.9215
s1 -1.8e+09 -80.9215
s1 -1.7e+09 -80.9215
s1 -1.6e+09 -80.9215
s1 -1.5e+09 -80.9215
s1 -1.4e+09 -80.9215
s1 -1.3e+09 -80.9215
s1 -1.2e+09 -80.9215
s1 -1.1e+09 -80.9215
s1 -1e+09 -80.9215
s1 -9e+08 -80.9215
s1 -8e+08 -80.9215
s1 -7e+08 -80.9215
s1 -6e+08 -80.9215
s1 -5e+08 -80.9215
s1 -4e+08 -80.9215
s1 -3e+08 -80.9215
s1 -2e+08 -80.9215
s1 -1e+08 -80.9215

Post by **ted** » Sun Jul 16, 2006 1:24 pm

The problem is that repeated executions do not yield the same results.
Eg, in the call sequence init(0.1), init(0.5), init(0.1), the first init(0.1) diverges while the second converges.

So something has changed between the first and second invocations of init(0.1). Most likely
the initial conditions are different because of an incorrect INITIAL block in one or more of the
many mechanisms, which results in one or more variables retaining values that are left over
from previous runs. Consequently the initial conditions of the entire system can vary from
run to run.

Another possibility is a change in the model specification itself. Are there any point
processes or sections that are attached to internal locations of another section? Your
init(0.5) discretizes the model with a coarser spatial grid, which will reposition any
internally-attached point process or section (except for those located at 0.5). When
you later execute init(0.1), those repositioned point processes and sections will not go
back to their original locations.

A final comment: init() happens to be the name of a key proc in NEURON's standard
run system (see chapter 8 of The NEURON Book). Redefining it in the way shown
here, which comingles it with code that specifies model parameters (the forall loop in
mulfit_cell_init() that specifies nseg and conductance densities) may make it very
difficult to use your model with the standard run system, complicating maintenance
and debugging of your code.

jaambros · Post by **jaambros** » Sun Jul 16, 2006 2:46 pm

I isolated the problem so it can be reproduced by others.
In the example below I use the pyramidal cell from the nrn demo and the hh mechanism.
Note that the first two calls to sinit(0.2) yield identical results.
However, after an intermediate nseg reset, sinit(0.2) yields different results.
Somehow altering nseg had side effects that were not finitialized().

This discrepancy makes it very difficult to compare spatial resolution effects in my project.

Code: Select all

load_file( "nrngui.hoc" )
// load pyramid.nrn from the Neuron demo
load_file( "/axon/d1/Users/jose/Software/src/Neuron/nrn/share/nrn/demo/pyramid.nrn" )

iss_idebug = 1
proc init(){ local dtsav, was_active, was_v
  if( iss_idebug ) printf( "iss t %g dt %g tstop %g cvode %g\n", t, dt, tstop, cvode.active() )
  finitialize( -70 )
  dtsav	= dt
  was_active  = cvode.active()
  t     = -1e9
  dt	= 1e8
  cvode.active( 0 )
  while( t < 0 ) { 
    if( iss_idebug ) printf( "s1 %g %g\n", t, v )
    fadvance()
  }
  if( iss_idebug ) printf( "iss Done s1 t %g v %g\n", t, v )

  cvode.active( was_active )		// restore cvode and dt state
  dt  = dtsav
  t   = 0
  if( cvode.active() ) {
    cvode.re_init()
  } else {
    fcurrent()
  }
  frecord_init()
}

proc set_nseg(){
  nseg_tot = 0
  soma area( 0.5 )
  forall { 
    nseg = int((L/($1 *lambda_f(100))+.9)/2)*2 + 1 
    nseg_tot += nseg
  }
  printf( "lambda-d %g nseg_tot %d\n", $1, nseg_tot )
}

proc mulfit_cell_init() { 
  forall { 
    Ra		= 31.3 
    cm		= 0.85
    g_pas	= 1 / ( 82 * 1000 )
  }
  set_nseg( $1 )
}

proc sinit() { 
  mulfit_cell_init( $1 )
  init( )
}

proc do_cell(){
  forall { 
    insert pas
    insert hh
    ion_style( "na_ion", 0, 2, 1, 0, 1 )
    ion_style( "k_ion", 0, 2, 1, 0, 1 )
  }
  nai0_na_ion = 10
  nao0_na_ion = 151
  ki0_k_ion = 140
  ko0_k_ion = 3.4
}

do_cell() 
  
tstop 		= 100
cvode.active(1)
cvode.atol(1.e-3)

access soma

sinit( 0.2 )
sinit( 0.2 )
forall nseg=1
sinit( 0.2 )

The output is

Code: Select all

NEURON -- Version 5.8 2005-10-7 13:46:29 Main (85)
by John W. Moore, Michael Hines, and Ted Carnevale
Duke and Yale University -- Copyright 1984-2005

nrniv: unable to open font "*Arial*bold*--12*", using "fixed"
	1 
	1 
	1 
	0.001 
lambda-d 0.2 nseg_tot 103
iss t 0 dt 0.025 tstop 100 cvode 1
s1 -1e+09 -70
s1 -9e+08 -63.3819
s1 -8e+08 -71.2783
s1 -7e+08 -61.8814
s1 -6e+08 -72.4312
s1 -5e+08 -60.5967
s1 -4e+08 -73.1389
s1 -3e+08 -59.893
s1 -2e+08 -73.4403
s1 -1e+08 -59.6181
iss Done s1 t 0 v -73.5265
lambda-d 0.2 nseg_tot 103
iss t 0 dt 0.025 tstop 100 cvode 1
s1 -1e+09 -70
s1 -9e+08 -63.3819
s1 -8e+08 -71.2783
s1 -7e+08 -61.8814
s1 -6e+08 -72.4312
s1 -5e+08 -60.5967
s1 -4e+08 -73.1389
s1 -3e+08 -59.893
s1 -2e+08 -73.4403
s1 -1e+08 -59.6181
iss Done s1 t 0 v -73.5265
lambda-d 0.2 nseg_tot 103
iss t 0 dt 0.025 tstop 100 cvode 1
s1 -1e+09 -70
s1 -9e+08 -64.3795
s1 -8e+08 -72.1959
s1 -7e+08 -61.4533
s1 -6e+08 -74.8538
s1 -5e+08 -58.7143
s1 -4e+08 -76.0981
s1 -3e+08 -57.7707
s1 -2e+08 -76.2847
s1 -1e+08 -57.6484
iss Done s1 t 0 v -76.2957
oc>

Post by **ted** » Mon Jul 17, 2006 11:08 am

The problem occurs even with a simple stylized (L, diam) ball and stick model, and
discretization by the d_lambda rule is irrelevant--just changing nseg to 1 and then back to
its original value will bring it out. The problem appears to be related to the initialization of
ena and ek when ion_style is called with the particular combination of arguments used
here--it disappears if the ion_style statements are commented out and the numerical
values of ena and ek are asserted by simple assignment.

Post by **hines** » Sun Jul 23, 2006 12:57 pm

I've fixed the error in which new segments do not inherit a previously specified ion_style but end up having the default ion style.
http://miris.med.yale.edu/cgi-bin/trac. ... geset/1462

Be aware that the only guaranteed lossless transformation is increasing nseg by an odd factor then decreasing to the original nseg. When nseg is increased, continuous inhomogeneities are not automatically recomputed, instead all range variables that are added get the values of the old segment they are contained in. This is fine for constants but one should consider re-executing the code that sets range variables after an nseg change if continuous functions of range variables are involved. The CellBuilder in "Continuous Create" mode will recompute the inhomogeneous range variables when dlambda is changed.However if one uses the hoc file derived from the CellBuilder and changes the segmentation then the proc biophys_inhomo() for that instance should be called.