Global structure, robustness, and modulation of neuronal models (Goldman et al. 2001)

 Download zip file 
Help downloading and running models
Accession:127878
"The electrical characteristics of many neurons are remarkably robust in the face of changing internal and external conditions. At the same time, neurons can be highly sensitive to neuromodulators. We find correlates of this dual robustness and sensitivity in a global analysis of the structure of a conductance-based model neuron. ..."
Reference:
1 . Goldman MS, Golowasch J, Marder E, Abbott LF (2001) Global structure, robustness, and modulation of neuronal models. J Neurosci 21:5229-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: Stomatogastric ganglion;
Cell Type(s):
Channel(s): I Na,t; I L high threshold; I T low threshold; I A; I K; I K,Ca; I Potassium;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: C or C++ program;
Model Concept(s): Bursting; Temporal Pattern Generation; Phase Response Curves;
Implementer(s): Goldman, Mark [msgoldman at ucdavis.edu];
Search NeuronDB for information about:  I Na,t; I L high threshold; I T low threshold; I A; I K; I K,Ca; I Potassium;
//random definitions I find useful, random number generators
#include <StdAfx.h>
#include <math.h>
#include <iostream>
#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include "mark_math.h"


//RANDOM NUMBER GENERATORS
//uniform dist. rand between 0.0 and 1.0(excluding endpts)
double ran1(long *idum)
{
	int j;
	long k;
	static long iy=0;
	static long iv[NTAB];
	double temp;

	if (*idum <= 0 || !iy) {		//Initialize
		if (-(*idum) < 1) *idum=1;  //be sure to prevent idum=0
		else *idum = -(*idum);
		for (j=NTAB+7; j>=0; j--) { //load the shuffle table (after 8 warm-ups)
			k=(*idum)/IQ;
			*idum=IA*(*idum-k*IQ)-IR*k;
			if (*idum < 0) *idum += IM;
			if (j < NTAB) iv[j] = *idum;
		}
		iy=iv[0];
	}
	k=(*idum)/IQ;							//start here when not initializing
	*idum=IA*(*idum-k*IQ)-IR*k;				//Compute idum=(IA*idum) % IM without overflows by Schrage's method
	if (*idum < 0) *idum += IM;
	j=iy/NDIV;								//will be in the range 0..NTAB-1
	iy=iv[j];								//output previously stored value and refill the shuffle table.
	iv[j] = *idum;
	if ((temp=AM*iy) > RNMX) return RNMX;	//Because users don't expect endpoint values
	else return temp;
}


//exp. dist. rand w/unit mean, using ran1(idum) as source of uniform deviates
double RandExp(long *idum) 
{
	double ran1(long *idum);
	double dum;

	do 
		dum=ran1(idum);
	while (dum == 0.0);
	return -log(dum);
}


//normal dist. rand w/mean zero, variance 1, using ran1(idum) 
double RandGauss(long *idum)
{
	
	//float ran1(long *idum);
	static int iset = 0;
	static double gset;
	double fac, rsq, v1, v2;

	if (iset == 0) {
		do {
			v1 = 2.0 * ran1(idum) - 1.0;
			v2 = 2.0 * ran1(idum) - 1.0;
			rsq = v1*v1 + v2*v2;
		}
		while (rsq >= 1.0 || rsq == 0.0);
		fac = sqrt(-2.0*log(rsq)/rsq);
		gset = v1*fac;
		iset = 1;
		return v2*fac;
	}
	else {
		iset = 0;
		return gset;
	}
}


//DATA TYPES -- matrices of doubles for number-crunching (from num. recipes)
void nrerror(char error_text[]) { //Numerical recipes standard error handler
	fprintf(stderr, "Numerical Recipes run-time error...\n");
	fprintf(stderr, "%s\n",error_text);
	fprintf(stderr, "...now exiting to system...\n");
	exit(1);
}

double *vector(long nl, long nh) {
//allocate a double vector with subscript range v[nl..nh]
	double *v;

	v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
	if (!v) nrerror("allocation failure in vector()");
	return v-nl+NR_END;
}

void free_vector(double *v, long nl, long nh) {
//free a double vector allocated with vector()
	free((FREE_ARG) (v+nl-NR_END));
}

double **matrix(long nrl, long nrh, long ncl, long nch)
// allocate a double matrix with subscript range m[nrl..nrh][ncl..nch]
{
	long i, nrow = nrh - nrl + 1, ncol = nch - ncl + 1;
	double **m;

	//allocate pointers to row
	m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
	if (!m) printf("allocation failure 1 in matrix()--whatever the hell that means\n");
	m += NR_END;
	m -= nrl;

	//allocate rows and set pointers to them 
	m[nrl] = (double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
	if (!m[nrl]) nrerror("allocation failure 2 in matrix()--whatever the hell that means");
	m[nrl] += NR_END;
	m[nrl] -= ncl;

	for (i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

	//return pointer ro array of pointers to rows
	return m;
}

void free_matrix(double **m, long nrl, long nrh, long ncl, long nch)
// free a double matrix allocated by matrix()
{
	free((FREE_ARG) (m[nrl] + ncl - NR_END));
	free((FREE_ARG) (m + nrl - NR_END));
}

//USEFUL FUNCTIONS
//returns max of a and b
double Max(double a, double b) {
	if (a > b)
		return(a);
	else return(b);
}

int Factorial(int n) {
	int nFactorial;
	if (n < 0) AfxMessageBox ("Factorial of number < 0 requested!");
	if (n < 2) { 
		nFactorial = 1;
	} else {
		nFactorial = 1;
		for (int i = 2; i <= n; i++) {
			nFactorial *= i;
		}
	}
	return(nFactorial);
}

double Poisson(int N, double AveN) { //value of Poisson distribution P(N)
	return(pow(AveN,N)*exp(-AveN)/Factorial(N));
}

double Binomial(int n, int N, double p) { //value of Binomial dist. P(n,N), n<N
	if (n > N) AfxMessageBox("Binomial n is > N");
	return(Factorial(N)*pow(p,n)*pow(1-p,(N-n))/(Factorial(n)*Factorial(N-n)));
}

//returns mean of elements of vector
double Mean(double *Vect, long NumElements) {
	if (NumElements < 1) {
		AfxMessageBox("0 elements: can't compute mean");
		return(-1.0);
	} else {
		double Ave = 0.0;
		for (long i = 0; i < NumElements; i++) {
			Ave += Vect[i];
		}
		Ave /= (double)NumElements;
		return(Ave);
	}
}

double StdDev(double *Vect, long NumElements) {
	if (NumElements < 2) return(-1.0);
	else {
		double Sigma;
		double Variance;
		double Sum = 0.0;
		double SquaresSum = 0.0;
		for (long i = 0; i < NumElements; i++) {
			Sum += Vect[i];
			SquaresSum += Vect[i]*Vect[i];
		}
		Variance = (SquaresSum - (Sum*Sum/(double)NumElements))/((double)NumElements - 1.0);
		Sigma = sqrt(Variance);
		return(Sigma);
	}
}

//Numerical recipes moment/cumulant computer: but I index from 0 to n-1, not 1 to n
void moment(double *data, long n, double *ave, double *adev, double *sdev, double *var, double *skew, double *curt) {
//Given an array of data, this routine return mean "ave", average deviation "adev", standard deviation "sdev"
//variance "var", skewness "skew", and kurtosis "curt"
	void nrerror(char error_text[]);
	long j;
	double ep = 0.0, s, p;

	if (n <= 1) nrerror("n must be at least 2 in moment");
	s = 0.0;								//1st pass to get the mean
	for (j=0; j < n; j++) s += data[j];
	*ave = s/n;
	*adev=(*var)=(*skew)=(*curt)=0.0;		//2nd pass to get the first (absolute), 2nd, 3rd, & 4th moments
	for (j=0; j < n; j++) {					//of the deviation from the mean
		*adev += fabs(s=data[j]-(*ave));
		ep += s;
		*var += (p=s*s);
		*skew += (p *= s);
		*curt += (p *= s);
	}
	*adev /= n;
	*var=(*var-ep*ep/n)/(n-1);				//Corrected 2-pass formula
	*sdev=sqrt(*var);						//Put the pieces together according to the conventional def.'s
	if (*var) {
		*skew /= (n*(*var)*(*sdev));
		*curt=(*curt)/(n*(*var)*(*var))-3.0;
	} else nrerror("No skew/kurtosis when variance = 0 (in moment)");
}

double ISIComputer(double *ISIList, double *SpikeList, long NumSpikes) {
	double MaxISI = 0.0;
	for (int i = 0; i < NumSpikes - 1; i++) {
		ISIList[i] = SpikeList[i+1] - SpikeList[i];
		if (ISIList[i] > MaxISI) 
			MaxISI = ISIList[i];
	}
	return(MaxISI);
}

void Histogram(double MaxEntry, int NumBins, long ListLength, double *List, double *hist) {
//MaxEntry is the largest entry in the bin	
	double BinSize = MaxEntry/NumBins;
	for (int bin = 0; bin < NumBins; bin++) {
		hist[bin] = 0.0;
	}
	for (int i=0; i < ListLength; i++) {
		bin	= (int)(List[i]/BinSize);
		hist[bin] += 1.0;
	}

}

//Next routine multiplies a matrix times a vector and places resulting vector in result
void MatrixMult(double* Result, double** Matrix, double* Vector, int NumRows, int NumCols) {
	for (int i = 0; i < NumRows; i++) {
		Result[i] = 0.0;
		for (int j = 0; j < NumCols; j++) {
			Result[i] += Matrix[i][j]*Vector[j];
		}
	}
}

//next routine is hardwired for a 2x2 matrix right now
void Eigenvals(double* Lambda, double** M) {
	Lambda[0] = (0.5*((M[0][0] + M[1][1]) + sqrt(SQR(M[0][0] - M[1][1]) + 4.*M[0][1]*M[1][0])));
	Lambda[1] = (0.5*((M[0][0] + M[1][1]) - sqrt(SQR(M[0][0] - M[1][1]) + 4.*M[0][1]*M[1][0])));
}

//next routine is hardwired for a 2x2 matrix right now
//actually returns the change of basis (FROM eigenbasis) matrix S, whose columns are the eigenvects
//Note: eigenvects NOT normalized
void Eigenvects(double** S, double** M) {
	double Lambda[2];  //hardwired at 2 for now
	Eigenvals(Lambda, M);
	S[0][0] = -M[0][1];	
	S[0][1] = -M[0][1];
	S[1][0] = M[0][0] - Lambda[0];
	S[1][1] = M[0][0] - Lambda[1];
}

//receives matrix M, returns eigenvals Lambda and eigenvects as columns of change of basis
//matrix (FROM eigenbasis) S
//Note: eigenvects NOT normalized
void Eigensystem(double** S, double* Lambda, double** M) {
	Lambda[0] = (0.5*((M[0][0] + M[1][1]) + sqrt(SQR(M[0][0] - M[1][1]) + 4.*M[0][1]*M[1][0])));
	Lambda[1] = (0.5*((M[0][0] + M[1][1]) - sqrt(SQR(M[0][0] - M[1][1]) + 4.*M[0][1]*M[1][0])));
	S[0][0] = -M[0][1];	
	S[0][1] = -M[0][1];
	S[1][0] = M[0][0] - Lambda[0];
	S[1][1] = M[0][0] - Lambda[1];
}

//next routine is hardwired for a 2x2 matrix right now
double Determinant(double** Matrix) {
	return ((Matrix[0][0] * Matrix[1][1]) - (Matrix[0][1] * Matrix[1][0]));
}

//next routine is hardwired for a 2x2 matrix right now; receives M, puts inverse in Inverse
void Inverse(double** Inverse, double** M) {
	double det = Determinant(M);
	Inverse[0][0] = M[1][1]/det;
	Inverse[0][1] = -M[0][1]/det;
	Inverse[1][0] = -M[1][0]/det;
	Inverse[1][1] = M[0][0]/det;
}

Loading data, please wait...