Increased computational accuracy in multi-compartmental cable models (Lindsay et al. 2005)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:129149
Compartmental models of dendrites are the most widely used tool for investigating their electrical behaviour. Traditional models assign a single potential to a compartment. This potential is associated with the membrane potential at the centre of the segment represented by the compartment. All input to that segment, independent of its location on the segment, is assumed to act at the centre of the segment with the potential of the compartment. By contrast, the compartmental model introduced in this article assigns a potential to each end of a segment, and takes into account the location of input to a segment on the model solution by partitioning the effect of this input between the axial currents at the proximal and distal boundaries of segments. For a given neuron, the new and traditional approaches to compartmental modelling use the same number of locations at which the membrane potential is to be determined, and lead to ordinary differential equations that are structurally identical. However, the solution achieved by the new approach gives an order of magnitude better accuracy and precision than that achieved by the latter in the presence of point process input.
Reference:
1 . Lindsay AE, Lindsay KA, Rosenberg JR (2005) Increased computational accuracy in multi-compartmental cable models by a novel approach for precise point process localization. J Comput Neurosci 19:21-38 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s):
Channel(s): I Na,t; I K;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON; C or C++ program;
Model Concept(s): Methods;
Implementer(s):
Search NeuronDB for information about:  I Na,t; I K;
/
LindsayEtAl2005
readme.txt
03-192.pdf
AnalyseResults.c
BitsAndPieces.c
CellData.dat
CompareSpikeTrain.c
Ed04.tex
ExactSolution.dat
GammaCode
Gen.tex
Gen1.tex
Gen2.tex
Gen3.tex
Gen4.tex
Gen5.tex
Gen6.tex
GenCom.c
GenCom1.c
GenCom2.c
GenComExactSoln.c
GenerateInput.c
GenerateInputText.c
GenRan.ran
GetNodeNumbers.c
Info100.dat
Info20.dat
Info200.dat
Info30.dat
Info300.dat
Info40.dat
Info400.dat
Info50.dat
Info500.dat
Info60.dat
Info70.dat
Info80.dat
Info90.dat
InputCurrents.dat
InputDendrite.dat
JaySpikeTrain.c
JayTest1.dat
JayTest100.dat
KenSpikeTrain.c
KenTest1.dat *
KenTest10.dat
KenTest100.dat *
KenTest10p.dat
KenTest1p.dat *
KenTest2.dat
KenTest2p.dat
KenTest3.dat
KenTest3p.dat
KenTest4.dat
KenTest4p.dat
KenTest5.dat
KenTest5p.dat
KenTest6.dat
KenTest6p.dat
KenTest7.dat
KenTest7p.dat
KenTest8.dat
KenTest8p.dat
KenTest9.dat
KenTest9p.dat
LU.c
Mean50.dat
Mean500.dat
mosinit.hoc
NC.pdf
NC.tex
NC1.tex
NC2.tex
NC3.tex
NC4.tex
NC5.tex
NC6.tex
NCFig2.eps *
NCFig3.eps *
NCFig4.eps *
NCFig5a.eps *
NCFig5b.eps *
NCFig6.eps *
NCPics.tex
NeuronDriver.hoc
NewComExactSoln.c
NewComp.pdf
NewComp.ps
NewComp.tex
NewComp.toc
NewComp1.tex
NewComp2.tex
NewComp3.tex
NewComp4.tex
NewComp5.tex
NewComp6.tex
NewCompFig1.eps
NewCompFig2.eps *
NewCompFig3.eps *
NewCompFig4.eps *
NewCompFig5a.eps *
NewCompFig5b.eps *
NewCompFig6.eps *
NewCompPics.tex
NewComSpikeTrain.c
NewRes.dat
NewRes60.dat
NewRes70.dat
NewRes80.dat
NewSynRes40.dat
NewTestCell.d3
NResults.res
OldComExactSoln.c
out.res
principles_01.tex
rand
Ratio.dat
RelErr.dat
ReviewOfSpines.pdf
SpikeTimes.dat
TestCell.d3
TestCell1.d3
TestCell2.d3
TestCell3.d3
TestCell4.d3
testcellnew2.hoc
TestCGS.c
TestGen1.c
TestSim.hoc
TestSim020.hoc
TestSim030.hoc
TestSim040.hoc
TestSim050.hoc
TestSim060.hoc
TestSim070.hoc
TestSim080.hoc
TestSim090.hoc
TestSim1.hoc
TestSim100.hoc
TestSim200.hoc
TestSim300.hoc
TestSim400.hoc
TestSim500
TestSim500.hoc
                            
FUNCTION ran_gamma(s, first) RESULT(fn_val)
USE random
IMPLICIT NONE

! Uses the algorithm in
! Marsaglia, G. and Tsang, W.W. (2000)

! Generates a random gamma deviate for shape parameter s >= 1.

REAL, INTENT(IN)    :: s
LOGICAL, INTENT(IN) :: first
REAL                :: fn_val

! Local variables
REAL, SAVE  :: c, d
REAL        :: u, v, x

IF (s < 1.0) THEN
   WRITE(*, *) 'Shape parameter must be >= 1'
   STOP
END IF

IF (first) THEN
  d = s - 1./3.
  c = 1.0/SQRT(9.0*d)
END IF

! Start of main loop
DO

! Generate v = (1+cx)^3 where x is random normal; repeat if v <= 0.

  DO
    x = random_normal()
    v = (1.0 + c*x)**3
    IF (v > 0.0) EXIT
  END DO

! Generate uniform variable U

  CALL RANDOM_NUMBER(u)
  IF (u < 1.0 - 0.0331*x**4) THEN
    fn_val = d*v
    EXIT
  ELSE IF (LOG(u) < 0.5*x**2 + d*(1.0 - v + LOG(v))) THEN
    fn_val = d*v
    EXIT
  END IF
END DO

RETURN
END FUNCTION ran_gamma



PROGRAM test_rgamma
USE random
IMPLICIT NONE

INTERFACE
  FUNCTION ran_gamma(s, first) RESULT(fn_val)
    USE random
    IMPLICIT NONE
    REAL, INTENT(IN)    :: s
    LOGICAL, INTENT(IN) :: first
    REAL                :: fn_val
  END FUNCTION ran_gamma
END INTERFACE

REAL, PARAMETER  :: shape(8) = (/ 1.0, 2.0, 3.0, 4.0, 5.0, 10.0, 20.0, 50.0 /)
INTEGER          :: i, j
REAL             :: aver, dev, sumsq, start, finish, x

DO i = 1, 8
  WRITE(*, *)
  WRITE(*, '(a, f7.1)') ' Shape parameter = ', shape(i)

! First time the new algorithm

  CALL CPU_TIME(start)
  x = ran_gamma(shape(i), .TRUE.)
  aver = x
  sumsq = 0.0
  DO j = 2, 1000000
    x = ran_gamma(shape(i), .FALSE.)
    dev = x - aver
    aver = aver + dev/j
    sumsq = sumsq + dev*(x-aver)
  END DO
  CALL CPU_TIME(finish)
  sumsq = SQRT(sumsq/999999.)
  WRITE(*, *) 'Using Marsaglia & Tsang algorithm'
  WRITE(*, '(a, f7.2, a, 2f7.3)')  &
        ' Time = ', finish - start, 'sec.   Mean & st.devn. = ', aver, sumsq

! Repeat using algorithm in module random

  CALL CPU_TIME(start)
  x = random_gamma(shape(i), .TRUE.)
  aver = x
  sumsq = 0.0
  DO j = 2, 1000000
    x = random_gamma(shape(i), .FALSE.)
    dev = x - aver
    aver = aver + dev/j
    sumsq = sumsq + dev*(x-aver)
  END DO
  CALL CPU_TIME(finish)
  sumsq = SQRT(sumsq/999999.)
  WRITE(*, *) 'Using algorithm in module random'
  WRITE(*, '(a, f7.2, a, 2f7.3)')  &
        ' Time = ', finish - start, 'sec.   Mean & st.devn. = ', aver, sumsq
END DO

STOP
END PROGRAM test_rgamma


DAGPUNAR PROG
-------------

FUNCTION random_gamma(s, b, first) RESULT(fn_val)

! Adapted from Fortran 77 code from the book:
!     Dagpunar, J. 'Principles of random variate generation'
!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

!     N.B. This version is in `double precision' and includes scaling

!     FUNCTION GENERATES A RANDOM GAMMA VARIATE.
!     CALLS EITHER random_gamma1 (S > 1.0)
!     OR random_exponential (S = 1.0)
!     OR random_gamma2 (S < 1.0).

!     S = SHAPE PARAMETER OF DISTRIBUTION (0 < REAL).
!     B = Scale parameter

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)

REAL (dp), INTENT(IN)  :: s, b
LOGICAL, INTENT(IN)    :: first
REAL (dp)              :: fn_val

! Local parameters
REAL (dp), PARAMETER  :: one = 1.0_dp, zero = 0.0_dp

IF (s <= zero) THEN
  WRITE(*, *) 'SHAPE PARAMETER VALUE MUST BE POSITIVE'
  STOP
END IF

IF (s >= one) THEN
  fn_val = random_gamma1(s, first)
ELSE IF (s < one) THEN
  fn_val = random_gamma2(s, first)
END IF

! Now scale the random variable
fn_val = b * fn_val
RETURN

CONTAINS


FUNCTION random_gamma1(s, first) RESULT(fn_val)

! Adapted from Fortran 77 code from the book:
!     Dagpunar, J. 'Principles of random variate generation'
!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM
! A GAMMA DISTRIBUTION WITH DENSITY PROPORTIONAL TO GAMMA**(S-1)*EXP(-GAMMA),
! BASED UPON BEST'S T DISTRIBUTION METHOD

!     S = SHAPE PARAMETER OF DISTRIBUTION
!          (1.0 < REAL)

REAL (dp), INTENT(IN)  :: s
LOGICAL, INTENT(IN)    :: first
REAL (dp)              :: fn_val

!     Local variables
REAL (dp)             :: d, r, g, f, x
REAL (dp), SAVE       :: b, h
REAL (dp), PARAMETER  :: sixty4 = 64.0_dp, three = 3.0_dp, pt75 = 0.75_dp,  &
                         two = 2.0_dp, half = 0.5_dp

IF (s <= one) THEN
  WRITE(*, *) 'IMPERMISSIBLE SHAPE PARAMETER VALUE'
  STOP
END IF

IF (first) THEN                        ! Initialization, if necessary
  b = s - one
  h = SQRT(three*s - pt75)
END IF

DO
  CALL RANDOM_NUMBER(r)
  g = r - r*r
  IF (g <= zero) CYCLE
  f = (r - half)*h/SQRT(g)
  x = b + f
  IF (x <= zero) CYCLE
  CALL RANDOM_NUMBER(r)
  d = sixty4*g*(r*g)**2
  IF (d <= zero) EXIT
  IF (d*x < x - two*f*f) EXIT
  IF (LOG(d) < two*(b*LOG(x/b) - f)) EXIT
END DO
fn_val = x

RETURN
END FUNCTION random_gamma1



FUNCTION random_gamma2(s, first) RESULT(fn_val)

! Adapted from Fortran 77 code from the book:
!     Dagpunar, J. 'Principles of random variate generation'
!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM
! A GAMMA DISTRIBUTION WITH DENSITY PROPORTIONAL TO
! GAMMA2**(S-1) * EXP(-GAMMA2),
! USING A SWITCHING METHOD.

!    S = SHAPE PARAMETER OF DISTRIBUTION
!          (REAL < 1.0)

REAL (dp), INTENT(IN)  :: s
LOGICAL, INTENT(IN)    :: first
REAL (dp)              :: fn_val

!     Local variables
REAL (dp)             :: r, x, w
REAL (dp), SAVE       :: a, p, c, uf, vr, d
REAL (dp), PARAMETER  :: vsmall = EPSILON(one)

IF (s <= zero .OR. s >= one) THEN
  WRITE(*, *) 'SHAPE PARAMETER VALUE OUTSIDE PERMITTED RANGE'
  STOP
END IF

IF (first) THEN                        ! Initialization, if necessary
  a = one - s
  p = a/(a + s*EXP(-a))
  IF (s < vsmall) THEN
    WRITE(*, *) 'SHAPE PARAMETER VALUE TOO SMALL'
    STOP
  END IF
  c = one/s
  uf = p*(vsmall/a)**s
  vr = one - vsmall
  d = a*LOG(a)
END IF

DO
  CALL RANDOM_NUMBER(r)
  IF (r >= vr) THEN
    CYCLE
  ELSE IF (r > p) THEN
    x = a - LOG((one - r)/(one - p))
    w = a*LOG(x)-d
  ELSE IF (r > uf) THEN
    x = a*(r/p)**c
    w = x
  ELSE
    fn_val = zero
    RETURN
  END IF

  CALL RANDOM_NUMBER(r)
  IF (one-r <= w .AND. r > zero) THEN
    IF (r*(w + one) >= one) CYCLE
    IF (-LOG(r) <= w) CYCLE
  END IF
  EXIT
END DO

fn_val = x
RETURN

END FUNCTION random_gamma2

END FUNCTION random_gamma



PROGRAM demo_gamma
! Demo program to generate a small smaple of random gamma's

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)

REAL (dp)  :: shape, scale, rg(100), avge, mean
LOGICAL    :: first
INTEGER    :: i, order(100)

INTERFACE
  FUNCTION random_gamma(s, b, first) RESULT(fn_val)
    IMPLICIT NONE
    INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)
    REAL (dp), INTENT(IN)  :: s, b
    LOGICAL, INTENT(IN)    :: first
    REAL (dp)              :: fn_val
  END FUNCTION random_gamma

  RECURSIVE SUBROUTINE quick_sort(list, order)
    IMPLICIT NONE
    INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)
    REAL (dp), DIMENSION (:), INTENT(IN OUT)  :: list
    INTEGER, DIMENSION (:), INTENT(OUT)      :: order
  END SUBROUTINE quick_sort
END INTERFACE

WRITE(*, '(a)', ADVANCE='NO') ' Enter the shape & scale parameter values: '
READ(*, *) shape, scale

first = .TRUE.
DO i = 1, 100
  rg(i) = random_gamma(shape, scale, first)
  first = .FALSE.
END DO

CALL quick_sort(rg, order)
WRITE(*, '(" ", 10f7.3)') rg

mean = shape * scale
avge = SUM( rg ) / 100._dp
WRITE(*, '(a, g12.5, a, g12.5)')  &
         ' Population mean = ', mean, '  Sample mean = ', avge

STOP
END PROGRAM demo_gamma



RECURSIVE SUBROUTINE quick_sort(list, order)

! Quick sort routine from:
! Brainerd, W.S., Goldberg, C.H. & Adams, J.C. (1990) "Programmer's Guide to
! Fortran 90", McGraw-Hill  ISBN 0-07-000248-7, pages 149-150.
! Modified by Alan Miller to include an associated integer array which gives
! the positions of the elements in the original order.

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)

REAL (dp), DIMENSION (:), INTENT(IN OUT)  :: list
INTEGER, DIMENSION (:), INTENT(OUT)       :: order

! Local variable
INTEGER :: i

DO i = 1, SIZE(list)
  order(i) = i
END DO

CALL quick_sort_1(1, SIZE(list))
RETURN


CONTAINS


RECURSIVE SUBROUTINE quick_sort_1(left_end, right_end)

INTEGER, INTENT(IN) :: left_end, right_end

!     Local variables
INTEGER             :: i, j, itemp
REAL (dp)           :: reference, temp
INTEGER, PARAMETER  :: max_simple_sort_size = 6

IF (right_end < left_end + max_simple_sort_size) THEN
  ! Use interchange sort for small lists
  CALL interchange_sort(left_end, right_end)

ELSE
  ! Use partition ("quick") sort
  reference = list((left_end + right_end)/2)
  i = left_end - 1
  j = right_end + 1

  DO
    ! Scan list from left end until element >= reference is found
    DO
      i = i + 1
      IF (list(i) >= reference) EXIT
    END DO
    ! Scan list from right end until element <= reference is found
    DO
      j = j - 1
      IF (list(j) <= reference) EXIT
    END DO


    IF (i < j) THEN
      ! Swap two out-of-order elements
      temp = list(i)
      list(i) = list(j)
      list(j) = temp
      itemp = order(i)
      order(i) = order(j)
      order(j) = itemp
    ELSE IF (i == j) THEN
      i = i + 1
      EXIT
    ELSE
      EXIT
    END IF
  END DO

  IF (left_end < j) CALL quick_sort_1(left_end, j)
  IF (i < right_end) CALL quick_sort_1(i, right_end)
END IF

RETURN
END SUBROUTINE quick_sort_1


SUBROUTINE interchange_sort(left_end, right_end)

INTEGER, INTENT(IN) :: left_end, right_end

!     Local variables
INTEGER             :: i, j, itemp
REAL (dp)           :: temp

DO i = left_end, right_end - 1
  DO j = i+1, right_end
    IF (list(i) > list(j)) THEN
      temp = list(i)
      list(i) = list(j)
      list(j) = temp
      itemp = order(i)
      order(i) = order(j)
      order(j) = itemp
    END IF
  END DO
END DO

RETURN
END SUBROUTINE interchange_sort

END SUBROUTINE quick_sort

Loading data, please wait...