Region-specific atrophy in dendrites (Narayanan, Narayan, Chattarji, 2005)

 Download zip file 
Help downloading and running models
Accession:147756
...in this study, we develop an algorithm that uses statistics from precise morphometric analyses to systematically remodel neuronal reconstructions. We use the distribution function of the ratio of two normal distributed random variables to specify the probabilities of remodeling along various regions of the dendritic arborization. We then use these probabilities to drive an iterative algorithm for manipulating the dendritic tree in a region-specific manner. As a test, we apply this framework to a well characterized example of dendritic remodeling: stress-induced dendritic atrophy in hippocampal CA3 pyramidal cells. We show that our pruning algorithm is capable of eliciting atrophy that matches biological data from rodent models of chronic stress. <br>
Reference:
1 . Narayanan R, Narayan A, Chattarji S (2005) A probabilistic framework for region-specific remodeling of dendrites in three-dimensional neuronal reconstructions. Neural Comput 17:75-96 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell; Dendrite;
Brain Region(s)/Organism: Hippocampus;
Cell Type(s): Hippocampus CA3 pyramidal cell;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: C or C++ program;
Model Concept(s): Influence of Dendritic Geometry;
Implementer(s): Narayanan, Rishikesh [rishi at iisc.ac.in];
Search NeuronDB for information about:  Hippocampus CA3 pyramidal cell;
#include"Prune.h"

/*****************************************************************/
// Polar Method (Box, Muller and Marsaglia) for two independent
// normal distributions
/*****************************************************************/

twodouble Prune :: normal (double mu1, double var1, double mu2,
										double var2)
{
	double v1, v2, s;
	do {
		v1 = 2 * drand48() - 1;
		v2 = 2 * drand48() - 1;
		s = v1*v1 + v2*v2;
	}
	while (s >= 1 || s < 1e-30);
	s = sqrt((-2*log(s))/s);
	v1 *= s;
	v2 *= s;
	twodouble x;
	x.x1=v1*sqrt(var1)+mu1;
	x.x2=v2*sqrt(var2)+mu2;
	return (x);
}


/*****************************************************************/
// Taylor Series expansion to find area under Normal distribution
// P(z\le x)=0.5 + {{1}\over{\sqrt{2\pi}}} \sum_{k=0}^\infty 
// {{(-1)^k x^{2*k+1}} \over {(2k+1) 2^k k!}}
// z is standard normal.
/*****************************************************************/

double Prune :: pnorm(double x) // Gaussian CDF
{
	if(x<-6.5) return 0;
	if(x>6.5) return 1;
	double factK=1;
	double sum=0.0;
	double term=1.0;
	int k=0;
	while(fabs(term)>exp(-23)){
		term=(1.0/sqrt(2.0*M_PI))*pow(-1,k)*pow(x,2*k+1)/
						((2*k+1)*pow(2,k)*factK);
		sum+=term;
		k++;
		factK*=k;
	}
	sum+=0.5;
	if(sum<1e-10) sum=0;
	return sum;
}

/*****************************************************************/
// Returns pdf of (a+x)/(b+y), where x and y are independent 
// standard normal distributed.
/*****************************************************************/

double Prune :: ratio(double a, double b, double t)
{
	double q,tf;
	q=(b+a*t)/sqrt(1+t*t);
	tf=(exp(-0.5*(a*a+b*b))/(M_PI*(1+t*t)))*(1+(q/dnorm(q))*
					(pnorm(q)-pnorm(0)));
	return tf;				
}

/*****************************************************************/

double Prune :: dnorm(double x) 	// Gaussian PDF
{
	return ((1.0/sqrt(2.0*M_PI))*exp(-x*x/2));
}

/*****************************************************************/
// The actual values of pruning are generated here
/*****************************************************************/


void Prune :: setZeta(char * basefilename)
{
	int i,k;
	struct timeval tv;
	struct timezone tz;
	gettimeofday(&tv,&tz);
	int SEED=tv.tv_usec;
	srand48(SEED);
	
	char * filename=new char[35]; 

	sprintf(filename,"%s.cb",basefilename);
    ifstream cbfile (filename);		// Control, BP
	sprintf(filename,"%s.cd",basefilename);
    ifstream cdfile (filename);		// Control, DL
	sprintf(filename,"%s.sb",basefilename);
    ifstream sbfile (filename);		// Stress, BP
	sprintf(filename,"%s.sd",basefilename);
    ifstream sdfile (filename);		// Stress, DL

    if(!cbfile || !cdfile || !sbfile || !sdfile){
        cerr << "\nOne or More of the Stats files not found.\n";
        exit(1);
	}	

// ************* CONVENTION *************/
// xyzp; x \in {c,s} - control and stress
//       y \in {b,d} - B.P. or D.L.
// 		 z \in {a,b} - apical or basal
// 		 p \in {m,v} - mean or variance
// *************************************/

	double * cdam; double * cdav; double * cdbm; double * cdbv;	
	double * cbam; double * cbav; double * cbbm; double * cbbv;

	double * sdam; double * sdav; double * sdbm; double * sdbv;
	double * sbam; double * sbav; double * sbbm; double * sbbv;

	cdam=new double[antra]; cdav=new double[antra];
	cbam=new double[antra]; cbav=new double[antra];
	cdbm=new double[bntra]; cdbv=new double[bntra];
	cbbm=new double[bntra]; cbbv=new double[bntra];
	sdam=new double[antra]; sdav=new double[antra];
	sbam=new double[antra]; sbav=new double[antra];
	sdbm=new double[bntra]; sdbv=new double[bntra];
	sbbm=new double[bntra]; sbbv=new double[bntra];

	for (i=0; i<antra; i++){	// Apical means
		cdfile >> cdam[i]; cbfile >> cbam[i];
		sdfile >> sdam[i]; sbfile >> sbam[i];
	}	
	for (i=0; i<antra; i++){	// Apical variances
		cdfile >> cdav[i]; cbfile >> cbav[i];
		sdfile >> sdav[i]; sbfile >> sbav[i];
	}	
	for (i=0; i<bntra; i++){	// Basal means
		cdfile >> cdbm[i]; cbfile >> cbbm[i];
		sdfile >> sdbm[i]; sbfile >> sbbm[i];
	}	
	for (i=0; i<bntra; i++){	// Basal variances
		cdfile >> cdbv[i]; cbfile >> cbbv[i];
		sdfile >> sdbv[i]; sbfile >> sbbv[i];
	}	
	cdfile.close(); cbfile.close();
	sdfile.close(); sbfile.close();

	double a,b,max;

	double bpsum;
	double dlcnt, bpcnt;
	dlsum=bpsum=0.0;
	dlcnt=bpcnt=0.0;

	abpzeta=new double[sdata.na];
	bbpzeta=new double[sdata.nb];
	adlzeta=new double[sdata.na];
	bdlzeta=new double[sdata.nb];

//  Find maximum ratio on Apical side for D.L. and B.P. //
//  Generate Zeta on Apical Side for DL and BP //

	for(i=0; i<sdata.na; i++){
		if(i>(antra-1)) k=antra-1;
		else k=i;
		max=maxratio(sdam[k], sdav[k], cdam[k], cdav[k]);
		adlzeta[i]=generateZeta(sdam[k], sdav[k], cdam[k], cdav[k], max);
		adlzeta[i] *= sdata.apicalcount[i];

		max=maxratio(sbam[k], sbav[k], cbam[k], cbav[k]);
		abpzeta[i]=generateZeta(sbam[k], sbav[k], cbam[k], cbav[k], max);
		abpzeta[i] *= BPStats.apicalcount[i];

		dlsum += adlzeta[i] ;
		bpsum += abpzeta[i];
		dlcnt += sdata.apicalcount[i];
		bpcnt += BPStats.apicalcount[i];

	}	
	

//  Find maximum ratio on Basal side for DL and BP  //
//  Generate Zeta on Basal Side for DL and BP //

	for(i=0; i<sdata.nb; i++){
		if(i>(bntra-1)) k=bntra-1;
		else k=i;
		max=maxratio(sdbm[k], sdbv[k], cdbm[k], cdbv[k]);
		bdlzeta[i]=generateZeta(sdbm[k], sdbv[k], cdbm[k], cdbv[k], max);
		bdlzeta[i] *= sdata.basalcount[i];

		max=maxratio(sbbm[k], sbbv[k], cbbm[k], cbbv[k]);
		bbpzeta[i]=generateZeta(sbbm[k], sbbv[k], cbbm[k], cbbv[k], max);
		bbpzeta[i] *= BPStats.basalcount[i];

		dlsum += bdlzeta[i] ;
		bpsum += bbpzeta[i];
		dlcnt += sdata.basalcount[i];
		bpcnt += BPStats.basalcount[i];

	}	

/*
	cerr << "\nReduction in DL: " << dlsum ;
	cerr << "\nTotal DL:        " << dlcnt ;
	cerr << "\nReduction in BP: " << bpsum ;
	cerr << "\nTotal BP:        " << bpcnt ;
	cerr << endl ;
*/

	delete(cdam); delete(cdav); 
	delete(cbam); delete(cbav); 
	delete(cdbm); delete(cdbv); 
	delete(cbbm); delete(cbbv); 
	delete(sdam); delete(sdav); 
	delete(sbam); delete(sbav); 
	delete(sdbm); delete(sdbv); 
	delete(sbbm); delete(sbbv); 
}

/*****************************************************************/

double Prune :: generateZeta(double sm, double sv, double cm, 
							double cv, double maxratio)
{							 
	if(!maxratio) return drand48();

	double vara=1,varb=1;
	double a,b;
	if (sv) a=sm/sv; 			// Avoid division by zero
	else if (!sm) return 0;		// Variance is zero and mean is zero
	else {a=sm; vara=0;}		// Variance is zero and mean is nonzero
		
	if (cv) b=cm/cv;
	else if (!cm) return 0; 
	else {b=cm; varb=0;}
		
	twodouble zeta;
	double t;
	double dzeta;
	double min;
	if(sm>cm) min=1.0;
	else min=sm/cm;
	do{
		if(!vara && !varb) {
			return 1-sm/cm;
		}
		else if(!vara){
			zeta=normal(sm,0,b,1);
			t=zeta.x1/zeta.x2;
			dzeta=t/cv;  
		}	
		else if(!varb){
			zeta=normal(a,1,cm,0);
			t=zeta.x1/zeta.x2;
			dzeta=t*sv;
		}	
		else {
			zeta=normal(a,1,b,1);	// This is ORIG.
			//zeta=normal(sm,sv,cm,cv);			// TEST
			t=zeta.x1/zeta.x2;
			dzeta=(t*sv)/(cv);
		}	
	}
	while(ratio(a,b,t)<(0.1*maxratio) || dzeta > 1.0 || dzeta < 0.0);
	//while(dzeta > 1.0 || dzeta < (min*0.88));
	//while(ratio(a,b,t)<(0.65*maxratio) || dzeta > 1.0 || dzeta < (min*0.85));
	return (1-dzeta);
}	

/*****************************************************************/

double Prune :: maxratio (double sm, double sv, double cm, double cv)
{
	double a,b;
	if (sv) a=sm/sv; 			// Avoid division by zero
	else if (!sm) return 0;		// Variance is zero and mean is zero
	else a=sm;					// Variance is zero and mean is nonzero

	if (cv) b=cm/cv;
	else if (!cm) return 0; 
	else b=cm;

	double max=0.0;
	double t;
	for(t=-3; t<5; t+=0.01)
		if(max<ratio(a,b,t)) 
			max=ratio(a,b,t);
	return max;		
}	

/*****************************************************************/

void Prune :: setProbabilities()
{
	setDLProbs();
	setBPProbs();
}

/*****************************************************************/

void Prune :: setBPProbs()
{
	int i;
	double P;

	// Apical

	for(i=0; i<sdata.na; i++){
		if(abpzeta[i]==0 || adlzeta[i]==0)
			abpP[i]=0; 
		else{
			if(abpprune[i]==0) abpprune[i]=1;
			if(adlprune[i]==0) adlprune[i]=RESLN;
			P=(adlzeta[i]*abpprune[i])/(abpzeta[i]*adlprune[i]);
			abpP[i]=pow(adlP[i],P);
		}
	}	

	// Basal

	for(i=0; i<sdata.nb; i++){
		if(bbpzeta[i]==0 || bdlzeta[i]==0)
			bbpP[i]=0; 
		else{
			if(bbpprune[i]==0) bbpprune[i]=1;
			if(bdlprune[i]==0) bdlprune[i]=RESLN;
			P=(bdlzeta[i]*bbpprune[i])/(bbpzeta[i]*bdlprune[i]);
			bbpP[i]=pow(bdlP[i],P);
		}
	}	
}

/*****************************************************************/

void Prune :: setDLProbs()
{
/* Set Probabilities for DL Pruning by just setting the sum of  *
 * probabilities to 1 by normalizing (zeta-prune) values.		*/

	double dlSum=0.0;

	maxdlP=0;

	int i;
	for(i=0; i<sdata.na; i++){
		adlP[i]=adlzeta[i]-adlprune[i];
		dlSum+=adlP[i];
	}	
	
	for(i=0; i<sdata.nb; i++){
		bdlP[i]=bdlzeta[i]-bdlprune[i];
		dlSum+=bdlP[i];
	}	

	for(i=0; i<sdata.na; i++){
		adlP[i]/=dlSum;
		if(maxdlP<adlP[i]) maxdlP=adlP[i];
		//cerr << "ADLP: " << i << "   "  << adlP[i] << endl;
	}	

	for(i=0; i<sdata.nb; i++){
		bdlP[i]/=dlSum;
		if(maxdlP<bdlP[i]) maxdlP=bdlP[i];
		//cerr << "BDLP: " << i << "   "  << bdlP[i] << endl;
	}	

	maxdlP += 0.1*maxdlP;
}
/*****************************************************************/

Loading data, please wait...