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;
// ZhengModelDoc.cpp : implementation of the CZhengModelDoc class
//

#include "stdafx.h"
#include "fstream"
using namespace std;
#include "iostream"
#include "stdio.h"
#include "ZhengModelHeaders.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/////////////////////////////////////////////////////////////////////////////
// CZhengModelDoc

IMPLEMENT_SERIAL(CZhengModelDoc, CDocument, _VERSION_NUMBER)

BEGIN_MESSAGE_MAP(CZhengModelDoc, CDocument)
	//{{AFX_MSG_MAP(CZhengModelDoc)
		// NOTE - the ClassWizard will add and remove mapping macros here.
		//    DO NOT EDIT what you see in these blocks of generated code!
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CZhengModelDoc construction/destruction

CZhengModelDoc::CZhengModelDoc()
{
	m_dtmin = _DT;
	m_dtmax = _DT;
	m_dt = _DT;
	for (int i = 0; i < _MAX_COMPARTS; i++) {
		m_IInject[i] = 0.0;
		for (int step = 0; step < _MAX_INJ_TIMES; step++) {
			m_MaxInj[i][step] = 0.0;
		}
	}
	m_pTheNeuron = new CNeuron;
	m_bEndRun = FALSE;
}

CZhengModelDoc::~CZhengModelDoc()
{
	
	delete m_pTheNeuron;
}

BOOL CZhengModelDoc::OnNewDocument()
{
	if (!CDocument::OnNewDocument())
		return FALSE;

	return TRUE;
}


/////////////////////////////////////////////////////////////////////////////
// CZhengModelDoc diagnostics

#ifdef _DEBUG
void CZhengModelDoc::AssertValid() const
{
	CDocument::AssertValid();
}

void CZhengModelDoc::Dump(CDumpContext& dc) const
{
	CDocument::Dump(dc);
}
#endif //_DEBUG

/////////////////////////////////////////////////////////////////////////////
// CZhengModelDoc commands

BOOL CZhengModelDoc::OnOpenDocument(LPCTSTR lpszPathName) 
{	
	if (!CDocument::OnOpenDocument(lpszPathName)) {
		return FALSE;
	}
	
	
	return TRUE;
}

void CZhengModelDoc::DeleteContents() 
{
	/*if (m_bOpenCase) {
	/*for (int i=0; i<_MAX_COMPARTS; i++) {
		for (int j=0; j<_MAX_CHANNELS; j++) {
			delete GetCompart(i)->m_ChanArray[j];
		}
		GetCompart(i)->m_ChanArray.RemoveAll();
	}

		delete m_pTheNeuron;
	}*/
	
	CDocument::DeleteContents();
}


BOOL CZhengModelDoc::PeekAndPump()
{
	MSG msg;

	while (::PeekMessage (&msg, NULL, 0, 0, PM_NOREMOVE)) {
		if (!AfxGetApp()->PumpMessage()) {
			::PostQuitMessage(0);
			return FALSE;
		}
	}
	return TRUE;
}

void CZhengModelDoc::RunEngine(double tmax, BOOL SaveFlag, BOOL DisplayFlag, double InjOnTime, double InjOnTime2, double InjOffTime)
{
	//int num_steps = (int)(tmax/dt);
	//int num_steps_per_screen = (int)(_TIME_PER_SCREEN/dt); //number of calculated steps/screenful
	//int graph_frequency = num_steps_per_screen/_POINTS_PER_SCREEN; //integer division(drops remainder)
	//int gmax_update_frequency = (int)(_TIME_PER_GMAX_UPDATE/dt); //number of steps between updates
	
	double ScaleFactorNa = 1.0;
	double ScaleFactorCaT = 1.0;
	double ScaleFactorCaP = 1.0;
	double ScaleFactorA = 1.0;
	double ScaleFactorKCa = 1.0;
	double ScaleFactorKd = 1.0;
	double* SpikeTimes;
	SpikeTimes = vector(0, _MAX_SPIKES - 1);
	double* ISIList;
	ISIList = vector(0, _MAX_SPIKES - 1);
	double* ISIHistogram; //histogram of ISI's
	ISIHistogram = vector(0, _NUM_BINS - 1);
	double MeanISI, SigmaISI, CVISI;
	int MeanISIbin;
	double MaxISI = 0.0;
	//double ave, adev, sdev, var, skew, kurt;
	double AveRate = 0.0;
	double BurstRate = -1.0;
	double ZeroSpikeBurstRate = 0.0;
	double SpikesPerBurst = -1.0;
	long NumISISpike = 0;
	long NumISIBurst = 0;
	double MeanISIBurst = 0.0;
	double MeanISISpike = 0.0;
	double num_steps;
	double DCSensorAveVal = 0.0;
	double SlowSensorAveVal = 0.0;
	double FastSensorAveVal = 0.0;
	long NumSpikes = 0;
	int NumSilent = 0; //Number of runs classified as silent cells
	int NumBurst = 0;
	int NumTonic = 0;
	time_t timet;
	long seed = -(long)time(&timet);
	//long seed = _SEED;
	BOOL SpikeFlag = TRUE;
	BOOL JustSpikedFlag = FALSE;
	BOOL BurstSearchFlag = FALSE;
	BOOL BurstFlag = FALSE;
	double DataTakingTime;
	double TotalBurstAmplitude;
	double AmplitudePerBurst;
	double Vmin, VAve;
	BOOL NewBurstFlag = TRUE; //for 0-spike burster identification

	char str[5000];
	ofstream SummaryFile("SummaryData.dat");
	if (!SummaryFile) cerr << "couldn't open file StateGrid.dat" << endl;
	//ofstream SensorFile("SensStateDiagram.dat");
	//if (!SensorFile) cerr << "couldn't open sensor state diagram data file" << endl;
	ofstream VandIFile("VandIFile.dat");
	if (!VandIFile) cerr << "couldn't open file VandIFile.dat" << endl;
	ofstream PRCFile("PRC1.05HzEPSP2.dat"); //Phase response curve file with phase out vs. phase in of IPSP
	if (!PRCFile) cerr << "couldn't open file PRC.dat" << endl;

	//FOR PHASE RESPONSE CURVE COMPUTATION
	BOOL PRCDone;
	BOOL ComputingPhaseResponse;
	BOOL TestPulseFired;
	double PhaseStartTime,PhaseIn,PhaseOut,Period;
	int NumPhasePoints = 1;

	for (double IPSPSize = 0.0; IPSPSize < 0.0055; IPSPSize += 0.1) {
		GetSyn(_Soma, _Inhib)->m_gStep = IPSPSize;
	for (int ThisPhasePoint = 0; ThisPhasePoint < NumPhasePoints; ThisPhasePoint++) {
	PRCDone = FALSE;
	ComputingPhaseResponse = FALSE;
	TestPulseFired = FALSE;
	ScaleFactorNa = 1.0;
	ScaleFactorCaT = 1.0;
	ScaleFactorCaP = 1.0;
	ScaleFactorA = 1.0;
	ScaleFactorKCa = 1.0;
	ScaleFactorKd = 1.0;
	double gMinNa = 0.0;
	double gMinCaT = 0.0;
	double gMinCaP = 0.0;
	double gMinA = 0.0;
	double gMinKCa = 0.0;
	double gMinKd = 0.0;
	//for 1SBurster analysis
	double MeanNa = 283.15;
	double MeanCaT = 3.45;
	double MeanCaP = 2.76;
	double MeanA = 26.24;
	double MeanKCa = 145.60;
	double MeanKd = 37.95;
	double SigmaNa = 240.67;
	double SigmaCaT = 1.18;
	double SigmaCaP = 0.94;
	double SigmaA = 20.11;
	double SigmaKCa = 88.57;
	double SigmaKd = 50.62;
	double gNa;
	double gCaT; 
	double gA;
	double gKCa;
	double gKd; 
	double EllipseSum;
	//For covariance analysis (again of 1SB)
	double Means[5];
	Means[0] = MeanNa; Means[1] = MeanCaT; Means[2] = MeanA; Means[3] = MeanKCa; Means[4] = MeanKd;
	double StdDevs[5];
	StdDevs[0] = 242.16; StdDevs[1] = 88.53; StdDevs[2] = 43.15; StdDevs[3] = 19.83; StdDevs[4] = 1.10;
	double gMaxVect[5];
	double DotProduct;
	double **EigenMatrix;
	EigenMatrix = matrix(0, 4, 0, 4); //holds eigenvectors in columns
	EigenMatrix[0][0] = 0.9936; EigenMatrix[0][1] = -0.0147; EigenMatrix[0][2] = 0.1114; EigenMatrix[0][3] = -0.0057; EigenMatrix[0][4] = 0.0009;
	EigenMatrix[1][0] = -0.0007; EigenMatrix[1][1] = 0.0009; EigenMatrix[1][2] = -0.0024; EigenMatrix[1][3] = -0.0188; EigenMatrix[1][4] = 0.9998;
	EigenMatrix[2][0] = 0.0031; EigenMatrix[2][1] = 0.0179; EigenMatrix[2][2] = -0.0757; EigenMatrix[2][3] = -0.9968; EigenMatrix[2][4] = -0.0189;
	EigenMatrix[3][0] = -0.0136; EigenMatrix[3][1] = -0.9997; EigenMatrix[3][2] = -0.0119; EigenMatrix[3][3] = -0.0171; EigenMatrix[3][4] = 0.0005;
	EigenMatrix[4][0] = -0.1117; EigenMatrix[4][1] = -0.0090; EigenMatrix[4][2] = 0.9908; EigenMatrix[4][3] = -0.0758; EigenMatrix[4][4] = 0.0009;


	//FOR MANY RUN RANDOM
	for (int ThisRun = 0; ThisRun < 1; ThisRun++) {
	
	//FOR GRID
	
	//for (double NaGmax = 100.0; NaGmax < 750.0; NaGmax += 100.0) {
	/*
	GetChan(_Soma, _I_Na)->m_gMax = 200.0;//NaGmax;
	for (double CaTGmax = .625; CaTGmax < 4.5; CaTGmax += .625) {
		GetChan(_Soma, _I_CaT)->m_gMax = CaTGmax;
		GetChan(_Soma, _I_CaP)->m_gMax = 0.8*CaTGmax;
		for (double AGmax = 9.375; AGmax < 70.0; AGmax += 9.375) {
			GetChan(_Soma, _I_A)->m_gMax = AGmax;
			for (double KCaGmax = 37.5; KCaGmax < 275.0; KCaGmax += 37.5) {
				GetChan(_Soma, _I_KCa)->m_gMax = KCaGmax;
				for (double KdGmax = 25.0; KdGmax < 180.0; KdGmax += 25.0) {
					GetChan(_Soma, _I_Kd)->m_gMax = KdGmax;
	*/
		//initialize local counter for average sensor values
		num_steps = 0.0;
		DCSensorAveVal = 0.0;
		SlowSensorAveVal = 0.0;
		FastSensorAveVal = 0.0;
		//sprintf(str, "Run number %d", ThisRun);
		//AfxMessageBox(str);
		for (int i = 0; i < _MAX_SPIKES; i++) {
			SpikeTimes[i] = 0.0;
			ISIList[i] = 0.0;
		}
		AveRate = 0.0;
		NumSpikes = 0;
		//COMPUTE, DISPLAY, AND SAVE INITIAL VALUES, AND SET UP OUTPUT STREAMS
		//initialize V
		double V0[_MAX_COMPARTS];
		for (i = 0; i < _MAX_COMPARTS; i++) {
			V0[i] = GetCompart(i)->m_V0;
			GetCompart(i)->m_V = V0[i];
		}
		//initialize Ca (and corresponding reversal E's)
		double Ca0[_MAX_COMPARTS];
		for (i = 0; i < _MAX_COMPARTS; i++) {
			Ca0[i] = GetCompart(i)->m_Ca0;
			GetCompart(i)->m_Ca = Ca0[i];
			//update Ca channel reversal potentials
			for (int chan = 0; chan < _MAX_CHANNELS; chan++) {
				if (GetChan(i, chan)->m_bCaChannel) {
					if (GetChan(i, chan)->m_bExists) {
						GetChan(i, chan)->m_Update_E(Ca0[i]);
					}
				}
			}
		}
		//initialize gmax's
		//FOR RANDOM POINT RUN USING MANYRUNINITIAL3 AS INITIAL DATA
		//or, for 1 spike bursters, using ManyRunMean1SB as initial data
		
		for (i = 0; i < _MAX_COMPARTS; i++) {
			for (int chan = 0; chan < _MAX_CHANNELS; chan++) {
				switch (chan) {
					//commented lines are for DON'T VARY Na, Kd
					case _I_Na:  GetChan(i, chan)->m_gMax = gMinNa + ScaleFactorNa * GetChan(i, chan)->m_gMax0; break;
					//case _I_Na:  GetChan(i, chan)->m_gMax = gMinNa + 1.0 * GetChan(i, chan)->m_gMax0; break;
					case _I_A:  GetChan(i, chan)->m_gMax = gMinA + ScaleFactorA * GetChan(i, chan)->m_gMax0; break;
					case _I_KCa:  GetChan(i, chan)->m_gMax = gMinKCa + ScaleFactorKCa * GetChan(i, chan)->m_gMax0; break;
					case _I_Kd:  GetChan(i, chan)->m_gMax = gMinKd + ScaleFactorKd * GetChan(i, chan)->m_gMax0; break;
					//case _I_Kd:  GetChan(i, chan)->m_gMax = gMinKd + 1.0 * GetChan(i, chan)->m_gMax0; break;
					case _I_CaT:  GetChan(i, chan)->m_gMax = gMinCaT + ScaleFactorCaT * GetChan(i, chan)->m_gMax0; break;
					case _I_CaP:  GetChan(i, chan)->m_gMax = 0.8 * GetChan(i, _I_CaT)->m_gMax; break;
					default:  GetChan(i, chan)->m_gMax = GetChan(i, chan)->m_gMax0;
				}
			}
		}
		
		//
		
		//Ready the output streams
		/*sprintf(str, "Spikes_%.3g_%.3g_%.3g_%.2g_%.2g.dat", GetChan(_Soma, _I_A)->m_gMax, GetChan(_Soma, _I_KCa)->m_gMax, GetChan(_Soma, _I_Kd)->m_gMax, GetRunPage()->m_MaxInjSoma, GetRunPage()->m_MaxInjSoma2);
		ofstream SpikeFile(str);
		//AfxMessageBox(str);
		if (!SpikeFile) cerr << "couldn't open Spikes data file" << endl;*/
		//Voltage file saving next 3 lines: 
		/*sprintf(str, "V_%.3g_%.3g_%.3g_%.2g_%.2g.dat", GetChan(_Soma, _I_A)->m_gMax, GetChan(_Soma, _I_KCa)->m_gMax, GetChan(_Soma, _I_Kd)->m_gMax, GetRunPage()->m_MaxInjSoma, GetRunPage()->m_MaxInjSoma2);
		ofstream Vfile(str);
		 (!Vfile) cerr << "couldn't open Voltage data file" << endl;*/

		//make a "% error sphere" of radius 1 and reject values that aren't within it -- for next time through the loop
		//do {
		//RANDOM POINT RUN USING MANYRUNINITIAL3 AS INITIAL DATA
		
			ScaleFactorNa = 2.0 * ran1(&seed);
			ScaleFactorCaT = 10.0 * ran1(&seed);
			ScaleFactorCaP = ScaleFactorCaT;
			ScaleFactorA = 5.0 * ran1(&seed);
			ScaleFactorKCa = 6.0 * ran1(&seed);
			ScaleFactorKd = 10.0 * ran1(&seed);
		
		//RANDOM POINT RUN to find covariance of ellipses: cuz rejection quick, just start w/ManyRunInitial3
		//and then enforce fitting in 5-D covariance ellipse
		/*
		do {	
			ScaleFactorNa = 2.0 * ran1(&seed);
			ScaleFactorCaT = 10.0 * ran1(&seed);
			ScaleFactorCaP = ScaleFactorCaT;
			ScaleFactorA = 5.0 * ran1(&seed);
			ScaleFactorKCa = 6.0 * ran1(&seed);
			ScaleFactorKd = 10.0 * ran1(&seed);
			
			gMaxVect[0] = ScaleFactorNa * 400.0;
			gMaxVect[1] = ScaleFactorCaT * 0.5;
			gMaxVect[2] = ScaleFactorA * 15.0;
			gMaxVect[3] = ScaleFactorKCa * 50.0;
			gMaxVect[4] = ScaleFactorKd * 20.0;

			//gMaxVect[0] = GetChan(_Soma, _I_Na)->m_gMax;
			//gMaxVect[1] = GetChan(_Soma, _I_CaT)->m_gMax;
			//gMaxVect[2] = GetChan(_Soma, _I_A)->m_gMax;
			//gMaxVect[3] = GetChan(_Soma, _I_KCa)->m_gMax;
			//gMaxVect[4] = GetChan(_Soma, _I_Kd)->m_gMax;

			EllipseSum = 0.0;
			for (int eigen = 0; eigen <= 4; eigen++) {
				DotProduct = 0;
				for (int chan = 0; chan <= 4; chan++) {
					DotProduct += EigenMatrix[chan][eigen]*(gMaxVect[chan] - Means[chan]);
				}
				EllipseSum += DotProduct*DotProduct/(StdDevs[eigen]*StdDevs[eigen]);
			}
		sprintf(str, "Na = %g, EllipseSum = %g", gMaxVect[0], EllipseSum);
		AfxMessageBox(str);
		} while (EllipseSum > 1.0);
		//sprintf(str, "EllipseSum = %g", EllipseSum);
		//AfxMessageBox(str);
		*/

		//RANDOM POINT RUN USING ManyRunMean1SB as initial data
		/*
		do {
			gMinNa = Max(0, MeanNa - SigmaNa);
			gMinCaT = Max(0, MeanCaT - SigmaCaT);
			gMinA = Max(0, MeanA - SigmaA);
			gMinKCa = Max(0, MeanKCa - SigmaKCa);
			gMinKd = Max(0, MeanKd - SigmaKd);

			ScaleFactorNa = (MeanNa + SigmaNa - gMinNa) * ran1(&seed) / MeanNa;
			ScaleFactorCaT = (MeanCaT + SigmaCaT - gMinCaT) * ran1(&seed) / MeanCaT;
			ScaleFactorA = (MeanA + SigmaA - gMinA) * ran1(&seed) / MeanA;
			ScaleFactorKCa = (MeanKCa + SigmaKCa - gMinKCa) * ran1(&seed) / MeanKCa;
			ScaleFactorKd = (MeanKd + SigmaKd - gMinKd) * ran1(&seed) / MeanKd;
			
			gNa = gMinNa + ScaleFactorNa * MeanNa;
			gCaT = gMinCaT + ScaleFactorCaT * MeanCaT;
			gA = gMinA + ScaleFactorA * MeanA;
			gKCa = gMinKCa + ScaleFactorKCa * MeanKCa;
			gKd = gMinKd + ScaleFactorKd * MeanKd;

			EllipseSum = 0.0;
			EllipseSum += (gNa-MeanNa)*(gNa-MeanNa)/(SigmaNa*SigmaNa);
			EllipseSum += (gCaT-MeanCaT)*(gCaT-MeanCaT)/(SigmaCaT*SigmaCaT);
			EllipseSum += (gA-MeanA)*(gA-MeanA)/(SigmaA*SigmaA);
			EllipseSum += (gKCa-MeanKCa)*(gKCa-MeanKCa)/(SigmaKCa*SigmaKCa);
			EllipseSum += (gKd-MeanKd)*(gKd-MeanKd)/(SigmaKd*SigmaKd);
			//sprintf(str, "EllipseSum = %g", EllipseSum);
			//AfxMessageBox(str);
		} while (EllipseSum > 1.0);
		*/
		//} while ((SQR(ScaleFactorA - 1.0) + SQR(ScaleFactorKCa - 1.0) + SQR(ScaleFactorKd - 1.0)) > 1.0);
			//} while (ScaleFactorNa < 0.063); //ignore Na values < 25
	
		//PUT THIS HERE IF WANT FIRST POINT RANDOM ALSO
		/*
		for (i = 0; i < _MAX_COMPARTS; i++) {
			for (int chan = 0; chan < _MAX_CHANNELS; chan++) {
				switch (chan) {
					//commented lines are for DON'T VARY Na, Kd
					case _I_Na:  GetChan(i, chan)->m_gMax = gMinNa + ScaleFactorNa * GetChan(i, chan)->m_gMax0; break;
					//case _I_Na:  GetChan(i, chan)->m_gMax = gMinNa + 1.0 * GetChan(i, chan)->m_gMax0; break;
					case _I_A:  GetChan(i, chan)->m_gMax = gMinA + ScaleFactorA * GetChan(i, chan)->m_gMax0; break;
					case _I_KCa:  GetChan(i, chan)->m_gMax = gMinKCa + ScaleFactorKCa * GetChan(i, chan)->m_gMax0; break;
					case _I_Kd:  GetChan(i, chan)->m_gMax = gMinKd + ScaleFactorKd * GetChan(i, chan)->m_gMax0; break;
					//case _I_Kd:  GetChan(i, chan)->m_gMax = gMinKd + 1.0 * GetChan(i, chan)->m_gMax0; break;
					case _I_CaT:  GetChan(i, chan)->m_gMax = gMinCaT + ScaleFactorCaT * GetChan(i, chan)->m_gMax0; break;
					case _I_CaP:  GetChan(i, chan)->m_gMax = 0.8 * GetChan(i, _I_CaT)->m_gMax; break;
					default:  GetChan(i, chan)->m_gMax = GetChan(i, chan)->m_gMax0;
				}
			}
		}
		*/

		//start channel activation/inactivation at steady-state
		for (i = 0; i < _MAX_COMPARTS; i++) {
			for (int chan = 0; chan < _MAX_CHANNELS; chan++) {
				if ((GetChan(i, chan)->m_p) != 0)
					GetChan(i, chan)->m_m = GetChan(i, chan)->m_m_inf(V0[i], Ca0[i]);
				if ((GetChan(i, chan)->m_q) != 0)
					GetChan(i, chan)->m_h = GetChan(i, chan)->m_h_inf(V0[i]);
			}
		}
		//start sensor activation/inactivation at steady-state
		for (i = 0; i < _MAX_COMPARTS; i++) {
			for (int s = 0; s < _MAX_SENSORS; s++) {
				if ((GetSensor(i, s)->m_p) != 0)
					GetSensor(i, s)->m_m = GetSensor(i, s)->m_m_inf(GetCompart(i)->m_ICa());
				if ((GetSensor(i, s)->m_q) != 0)
					GetSensor(i, s)->m_h = GetSensor(i, s)->m_h_inf(GetCompart(i)->m_ICa());
			}
		}


		double t = 0.0;
		for (i = 0; i < _MAX_COMPARTS; i++) {
			m_IInject[i] = 0.0;
			if ((t >= InjOnTime) AND (t < InjOnTime2)) {
				m_IInject[i] = m_MaxInj[i][0]; //0=first step
				//NEXT LINES ADDED FOR CASE OF ADDING NOISE OF VARIANCE -IInject * m_C
				time_t timet;
				long seed = -(long)time(&timet);
				m_IInject[i] = (m_IInject[i]/sqrt(m_dt)) * RandGauss(&seed);
			}
			else if ((t >= InjOnTime2) AND (t < InjOffTime)) {
				m_IInject[i] = m_MaxInj[i][1]; //1=2nd step
				//NEXT LINES ADDED FOR CASE OF ADDING NOISE OF VARIANCE -IInject * m_C
				time_t timet;
				long seed = -(long)time(&timet);
				m_IInject[i] = (m_IInject[i]/sqrt(m_dt)) * RandGauss(&seed);
			}
		}
		if (DisplayFlag) {
			//GetPlotsDialog()->DoPlots((int)t, num_steps_per_screen, graph_frequency); //int(t) = t_step
			GetPlotsDialog()->DoPlots(t);
		}
		/*if (SaveFlag) {
			Vfile << t << " " << GetCompart(_Soma)->m_V << endl;//writing to file stuff for init vals
		}*/

		//RUN NEURON

		//ADAPTIVE STEP ALGORITHM--VARIABLE dt
		double VOld[_MAX_COMPARTS];
		double VNew[_MAX_COMPARTS];
		double num_gmax_updates = 0;
		BOOL step_up_flag;
		//if (dtmin != dtmax) {
			for (i = 0; i < _MAX_COMPARTS; i++) {
				VOld[i] = GetCompart(i)->m_V;
			}
			m_dt = m_dtmin;  //start at min for safety
			DataTakingTime = 20000.0;
			TotalBurstAmplitude = 0.0;
			ZeroSpikeBurstRate = 0.0;
			Vmin = 100.0;
			VAve = 0.0;
			while (t < tmax) {
				t += m_dt;
				//compute
				for (i = 0; i < _MAX_COMPARTS; i++) {
					m_IInject[i] = 0.0;
					if ((t >= InjOnTime) AND (t < InjOnTime2)) {
						m_IInject[i] = m_MaxInj[i][0]; //0=first step
						//NEXT LINES ADDED FOR CASE OF ADDING NOISE OF VARIANCE -IInject * m_C
						time_t timet;
						long seed = -(long)time(&timet);
						m_IInject[i] = (m_IInject[i]/sqrt(m_dt)) * RandGauss(&seed);
					}
					else if ((t >= InjOnTime2) AND (t < InjOffTime)) {
						m_IInject[i] = m_MaxInj[i][1]; //1=2nd step
						//NEXT LINES ADDED FOR CASE OF ADDING NOISE OF VARIANCE -IInject * m_C
						time_t timet;
						long seed = -(long)time(&timet);
						m_IInject[i] = (m_IInject[i]/sqrt(m_dt)) * RandGauss(&seed);
					}
				}		
				m_pTheNeuron->ComputeMe(m_IInject, m_dt);
				for (i=0; i < _MAX_COMPARTS; i++) {
					VNew[i] = GetCompart(i)->m_V;
				}
				
				//update step size
				for (i=0; i < _MAX_COMPARTS; i++) {
					if (fabs(VNew[i] - VOld[i]) > 1.0) { //decrease dt if any V's change by more than 1 mV
						m_dt = 0.2 * m_dt;
						if (m_dt < m_dtmin) {
							m_dt = m_dtmin;
						}
						break;
					}
				}
				step_up_flag = TRUE;
				for (i=0; i < _MAX_COMPARTS; i++) {
					if (fabs(VNew[i] - VOld[i]) > 0.5) { //increase dt only if ALL V's change by <= 0.2 mV
						step_up_flag = FALSE;	
					}
				}
				if (step_up_flag) {
					m_dt = 3.0 * m_dt;
					if (m_dt > m_dtmax) {
						m_dt = m_dtmax;
					}
				}

				//output
				if (DisplayFlag) {//for plotting fixed time & # of points/screen
					GetPlotsDialog()->DoPlots(t);
				}
				if ((int)(t/_TIME_PER_GMAX_UPDATE) != num_gmax_updates) { //currently every 500 ms
					num_gmax_updates = (int)(t/_TIME_PER_GMAX_UPDATE);
					GetNeuronPage()->UpdateGMaxDisplay();
					GetNeuronPage()->UpdateSensorValueDisplay();
				}
				//if (SaveFlag) {  //Hardwired to taking 20 sec's of data
					if ((t > tmax - DataTakingTime - 6000.0) AND (t < tmax - DataTakingTime - 1000.0)) {
						//search for minimum voltage, for case of zero-spike bursters
						if (VNew[_Soma] < Vmin) {
							Vmin = VNew[_Soma];
						}
					}
					if ((t > tmax - DataTakingTime - 1000.0) AND (t < tmax - 1000.0)) {
						num_steps += 1.0;
						//identify 0 spike bursters (throw this out if AveRate > 0.09)
						if ((NewBurstFlag) AND (VNew[_Soma] > _TRANSMIT_THRESHOLD) ) {
							NewBurstFlag = FALSE;
							//JustBurstFlag = TRUE;
							//NumBursts++;
							//will get overwritten if a spiking cell
							ZeroSpikeBurstRate += 1000.0/DataTakingTime; 
						}
						/*
						if ((NewBurstFlag) AND (VOld[_Soma] > VNew[_Soma])) {
							JustBurstFlag = FALSE;
							//BurstSearchFlag = TRUE;
						}*/
						if (((VNew[_Soma] - Vmin) < 2.0) AND (VNew[_Soma] < VOld[_Soma])) {
							NewBurstFlag = TRUE;
						}
						if ((SpikeFlag) AND (VNew[_Soma] > _SPIKE_THRESHOLD) AND (VOld[_Soma] < VNew[_Soma])) {
							SpikeFlag = FALSE;
							/*if (SaveFlag) {
								SpikeFile << t << endl;
							}*/
							SpikeTimes[NumSpikes] = t;

							//FOR COMPUTING PHASE RESPONSE CURVE--MAY REQUIRE CONSISTENT FIXED TIME STEP
							//HARDWIRED FOR 3 SPIKES/BURST AND FOR SPIKES W/IN BURST <100ms apart, BETWEEN > 100ms apart
							if (!PRCDone) {
								//possibly change below to several cycles of counting???
								if ((ComputingPhaseResponse) AND ((SpikeTimes[NumSpikes] - SpikeTimes[NumSpikes-1]) > 100.0)) {
									PRCDone = TRUE;
									PhaseOut = (SpikeTimes[NumSpikes] - PhaseStartTime - Period)/Period;
								}
								if ((SpikeTimes[NumSpikes] > 34000) AND (!ComputingPhaseResponse)
									AND ((SpikeTimes[NumSpikes] - SpikeTimes[NumSpikes - 1]) > 100.0)) {//start PRC calculation
									ComputingPhaseResponse = TRUE;
									Period = SpikeTimes[NumSpikes] - SpikeTimes[NumSpikes - 3];
									PhaseStartTime = SpikeTimes[NumSpikes];
								}
							}
							//

							JustSpikedFlag = TRUE;
							NumSpikes++;
							AveRate += 1000.0/DataTakingTime; //# of spikes/20 sec of data-taking
						}
						if ((JustSpikedFlag) AND (VOld[_Soma] > VNew[_Soma])) {
							JustSpikedFlag = FALSE;
							//BurstSearchFlag = TRUE;
						}
						if (VNew[_Soma] < _SPIKE_THRESHOLD) {
							SpikeFlag = TRUE;
						}
						//assume graded transmission for _trans_thresh < V, with saturation at trans_max
						if (VNew[_Soma] > _TRANSMIT_THRESHOLD) {
							TotalBurstAmplitude += (Max(VNew[_Soma],_TRANSMIT_MAX) - _TRANSMIT_THRESHOLD)*m_dt;
						}
						VAve += VNew[_Soma]*m_dt/DataTakingTime;
						/*
						if (VNew[_Soma] < _BURST_THRESHOLD) {//maybe redefine burst_thresh relative to min V?
							BurstSearchFlag = FALSE;
						}
						*/
						/*
						if ((BurstSearchFlag) AND (VOld[_Soma] < VNew[_Soma])) {
							BurstFlag = TRUE;
						}*/
					}
					//FOR PHASE RESPONSE CURVE COMPUTATION
					if (ComputingPhaseResponse) {
						if (!TestPulseFired) {
							//NOTE:NumSpikes-1 holds most recent spike time
							if ((t - SpikeTimes[NumSpikes-1]) > ((double)ThisPhasePoint*Period/(double)NumPhasePoints)) {
								GetSyn(_Soma, _Inhib)->m_SpikeReceived();
								PhaseIn = (t - PhaseStartTime)/Period;
								TestPulseFired = TRUE;
								//AfxMessageBox("Test pulse fired");
							}
						}
					}
					//
					if ((t > tmax - 11000.0) AND (t < tmax - 1000.0)) {
						//NEXT LINES ARE FOR 1-COMPARTMENT MAP OF SENSOR SPACE
						DCSensorAveVal += GetSensor(_Soma, _DC)->m_AveValue;
						SlowSensorAveVal += GetSensor(_Soma, _Slow)->m_AveValue;
						FastSensorAveVal += GetSensor(_Soma, _Fast)->m_AveValue;
					}
					//if ((t > InjOnTime - 21000.0) AND (t < InjOffTime - 16000.0)) {
					if ((t > InjOffTime - 19000.0) AND (t < InjOffTime - 16000.0)) {
						if (SaveFlag) {//voltages
							//PUT BACK IN NEXT LINE TO SAVE A VOLTAGE TRACE
							VandIFile << t << " " << VNew[_Soma] << endl; 
							/*
							VandIFile << " " << GetCompart(_Soma)->m_ITotal(0.0);
							for (int chan = _I_Leak; chan <= _I_h; chan++) {
								VandIFile << " " << GetChan(_Soma, chan)->m_I(VNew[_Soma]);
							}
							VandIFile << " " << GetCompart(_Soma)->m_Ca << endl;
							*/
						}
					}
				//}
				//re-set VOld
				for (i=0; i < _MAX_COMPARTS; i++) {
					VOld[i] = VNew[i];
				}
				if (!PeekAndPump())	// check for windows messages from user
					break;
				if (m_bEndRun) { //end run due to user clicking Finish button
					//m_bEndRun = FALSE;
					break;
				}
			}
			if (m_bEndRun) { //end run due to user clicking Finish button
				m_bEndRun = FALSE;
				break;
			}
			DCSensorAveVal /= num_steps;
			SlowSensorAveVal /= num_steps;
			FastSensorAveVal /= num_steps;
			MaxISI = ISIComputer(ISIList, SpikeTimes, NumSpikes);
			if (NumSpikes > 1) {
				Histogram(MaxISI, _NUM_BINS, NumSpikes - 1, ISIList, ISIHistogram);
				MeanISI = Mean(ISIList, NumSpikes - 1);
				SigmaISI = StdDev(ISIList, NumSpikes - 1);
				//next line requires NumSpikes >2
				//moment(ISIList, NumSpikes - 1, &ave, &adev, &sdev, &var, &skew, &kurt);
				CVISI = SigmaISI/MeanISI;
				sprintf(str, "Histogram: ");
				MeanISIbin = (int)((MeanISI/MaxISI)*(double)_NUM_BINS);
				for (i = 0; i < _NUM_BINS; i++) {
					sprintf(str, "%s %g", str, ISIHistogram[i]);
					if (i == MeanISIbin)
						sprintf(str, "%sR", str);
				}
				if (DisplayFlag) {
					AfxMessageBox(str);
				}
			} else {
				MeanISI = -1;
				SigmaISI = -1;
				CVISI = -1;
				//ave = -1;
				//sdev = -1;
				//kurt = -1;
			}
			//sprintf(str, "ISIs_%.3g_%.3g_%.3g_%.2g_%.2g.dat", GetChan(_Soma, _I_A)->m_gMax, GetChan(_Soma, _I_KCa)->m_gMax, GetChan(_Soma, _I_Kd)->m_gMax, GetRunPage()->m_MaxInjSoma, GetRunPage()->m_MaxInjSoma2);
			//ofstream ISIFile(str);
			//if (!ISIFile) cerr << "couldn't open ISI data file" << endl;
			if (SaveFlag) {
				//for (int i = 0; i < NumSpikes - 1; i++) {
					//ISIFile << ISIList[i] << endl;
				//}
				for (int chan = 0; chan < _MAX_CHANNELS; chan++) {
					SummaryFile << GetChan(_Soma, chan)->m_gMax << " ";
				}
				//for (sens = 0; sens < _MAX_SENSORS; sens++) {
					SummaryFile << DCSensorAveVal << " " << SlowSensorAveVal << " " << FastSensorAveVal << " "; 
				//}
			}
			if (AveRate < 0.09) {
				if (SaveFlag) {
					SummaryFile << "D "; //dead cell
				}
				if (ZeroSpikeBurstRate > 0.01) {
					AmplitudePerBurst = 1000.0*TotalBurstAmplitude/(ZeroSpikeBurstRate*DataTakingTime);
					BurstRate = ZeroSpikeBurstRate;
				} else {
					AmplitudePerBurst = -1.0;
					BurstRate = 0.0;
				}
				SpikesPerBurst = 0.0;
				NumSilent += 1;
			} else if (CVISI > 0.3) {//maybe better to split data in half about mean and use skewness
				if (SaveFlag) {
					SummaryFile << "B "; //bursting cell
				}
				NumBurst += 1;
				//calculate burst rate and #spikes/burst: crude, use average values of interburst and interspike ISI's
				NumISISpike = 0;
				NumISIBurst = 0;
				MeanISISpike = 0.0;
				MeanISIBurst = 0.0;
				for (i = 0; i < NumSpikes - 1; i++) {
					if (ISIList[i] <= MeanISI) {
						NumISISpike += 1;
						MeanISISpike += ISIList[i];
					} else if (ISIList[i] > MeanISI) {
						NumISIBurst += 1;
						MeanISIBurst += ISIList[i];
					} else AfxMessageBox("Error in calculating MeanISI");
				}
				MeanISISpike /= (double)(NumISISpike)*1000.0; //convert to [sec]
				MeanISIBurst /= (double)(NumISIBurst)*1000.0;
				BurstRate = (1.0 - AveRate*MeanISISpike)/(MeanISIBurst - MeanISISpike);
				AmplitudePerBurst = 1000.0*TotalBurstAmplitude/(BurstRate*DataTakingTime);
				SpikesPerBurst = AveRate/BurstRate;
				if (DisplayFlag) {
					sprintf(str, "AveRate %g, intraburst interspike %g, interburst %g, TransmitPerBurst %g \n%g Hz burster with %g spikes per burst", 
						    AveRate, MeanISISpike, MeanISIBurst, AmplitudePerBurst, BurstRate, SpikesPerBurst);
					AfxMessageBox(str);
					sprintf(str, "DC sensor %g, SlowSensor %g, FastSensor %g", DCSensorAveVal, SlowSensorAveVal, FastSensorAveVal);
					AfxMessageBox(str);
				}
			} else if (CVISI <= 0.3) {
				/*if (BurstFlag) {//single-spike burster
					if (DisplayFlag) {
						AfxMessageBox("1 spike burster");
					}
					if (SaveFlag) {
						SummaryFile << "B "; //bursting
					}
					BurstRate = AveRate;
					SpikesPerBurst = 1.0;
					NumBurst += 1;
				} else {//tonic	*/
					if (SaveFlag) {
						SummaryFile << "T "; //tonically firing
					}
					BurstRate = AveRate;
					AmplitudePerBurst = 1000.0*TotalBurstAmplitude/(AveRate*DataTakingTime);
					SpikesPerBurst = 1.0;
					NumTonic += 1;
					if (DisplayFlag) {
						sprintf(str, "AveRate %g, TransmitPerSpike %g", AveRate, AmplitudePerBurst);
						AfxMessageBox(str);
					}
				//}
			} else AfxMessageBox("Error in classifying dead, tonic, bursting");
			//Re-set flags
			JustSpikedFlag = FALSE;
			//JustBurstFlag = FALSE;
			BurstSearchFlag = FALSE;
			BurstFlag = FALSE;
			NewBurstFlag = FALSE;
			if (SaveFlag) {
				SummaryFile << AveRate << " " << BurstRate << " " << SpikesPerBurst << " ";
				//SummaryFile << ave << " " << sdev << " " << CVISI << " " << kurt << " " << seed << endl;
				SummaryFile << AmplitudePerBurst << " " << VAve << " " << CVISI << " " << seed << endl;
				//PRCFile << GetSyn(_Soma, _Inhib)->m_E << " " << GetSyn(_Soma, _Inhib)->m_tau_decay << " ";
				//PRCFile << GetSyn(_Soma, _Inhib)->m_gStep << " " << PhaseIn << " " << PhaseOut << endl;
			}
	//NEXT LINE FOR MANY RUN INITIAL
		} //for ThisRun
		} //for ThisPhasePoint
		} //for IPSPSize
	
   //FOR GRID
   /*
					} //for KdGmax
			} //for KCaGmax
		} //for AGmax
	} //for CaTGmax
	//} for NaGmax
   */
	if (DisplayFlag) {
		sprintf(str, "NumSilent = %d, NumBurst = %d, NumTonic = %d", NumSilent, NumBurst, NumTonic);
		AfxMessageBox(str);
	}
	//}  //end if (dtmin != dtmax)


	//free_vector(SpikeTimes, 0, _MAX_SPIKES - 1);
	//free_vector(ISIList, 0, _MAX_SPIKES - 1);
	//free_vector(ISIHistogram, 0, _NUM_BINS - 1);
	//free_matrix(EigenMatrix, 0, 4, 0, 4);
	
	//FOR FIXED TIME STEP dt
	/*int T;
	for (int t_step = 1; t_step <= num_steps; t_step++) {
		t = t_step * dt;
		T = t_step%num_steps_per_screen;
		for (i = 0; i < _MAX_COMPARTS; i++) {
			m_IInject[i] = 0.0;
			if ((t >= InjOnTime) AND (t < InjOffTime))
					m_IInject[i] = m_MaxInj[i];
		}
		m_pTheNeuron->ComputeMe(m_IInject, dt);
		if (DisplayFlag AND T%graph_frequency == 0) {//for plotting fixed time & # of points/screen
			GetPlotsDialog()->DoPlots(t_step, num_steps_per_screen, graph_frequency);
		}
		if (t_step%gmax_update_frequency == 0) {
			GetNeuronPage()->UpdateGMaxDisplay();
			GetNeuronPage()->UpdateSensorValueDisplay();
		}
		if (SaveFlag) {
			if (t_step%graph_frequency == 0)
				Vfile << t << " " << GetCompart(_Soma)->m_V << endl;//writing to file stuff for init vals
		}
		if (!PeekAndPump())	// check for windows messages from user
			break;
		if (m_bEndRun) { //end run due to user clicking Finish button
			m_bEndRun = FALSE;
			break;
		}
	}*/
}

/////////////////////////////////////////////////////////////////////////////
// CZhengModelDoc serialization

void CZhengModelDoc::Serialize(CArchive& ar)
{
	if (ar.IsStoring())
	{
		ar << m_dtmin;
		ar << m_dtmax;
		ar << m_dt;
		for (int i = 0; i < _MAX_COMPARTS; i++) {
			ar << m_IInject[i];
			for (int step = 0; step < _MAX_INJ_TIMES; step++) {
				ar << m_MaxInj[i][step];
			}
		}
		ar << m_bEndRun;
	}
	else
	{
		ar >> m_dtmin;
		ar >> m_dtmax;
		ar >> m_dt;
		for (int i = 0; i < _MAX_COMPARTS; i++) {
			ar >> m_IInject[i];
			for (int step = 0; step < _MAX_INJ_TIMES; step++) {
				ar >> m_MaxInj[i][step];
			}
		}
		ar >> m_bEndRun;

	}
	m_pTheNeuron->Serialize(ar);
	GetPlotsDialog()->Serialize(ar);
	GetRunPage()->Serialize(ar);
	GetNeuronPage()->Serialize(ar);
	GetDisplayPage()->Serialize(ar);
}

Loading data, please wait...