#include <iostream>
#include "nr3.h"
#include "odeint.h"
#include "ludcmp.h"
#include "stepper.h"
#include "stepperross.h"
#include <fstream>
#include <cmath>
using namespace std;
double a = 0.3646; // parameters for equation (1) in paper
double b = 8.6971;
double k1 = 2.0732;
double c = 0.4467;
double alpha = 10.0; // parameters for equation (2) in paper
double beta = 10.0;
double k2 = 0.2041;
double S = 0.1;
double gamma = 0.77; // "rise factor" (nM)
double eta = 0.91; // "decay rate" (/s)
double t_ox = 60.0; // "time constant for orexin" (s)
double t_ser = 10.0; // "time constant for serotonin/5_HT" (s)
double phi = 33.57; //"Release per stimulus rate" (nM)
double Km = 170.0; // (nM) Km Vmax Michaelis-Menton parameters
double Vmax = 1300; // (nM/s)
//NB Input/output (3) and (4) were substiuted directly subsitituted into the equations (1) and (2) in the paper
//dydx[0] = (-1/t_ox)*y[0] + (a/t_ox) + (b/t_ox)*(1/(1 + exp(k1/c)*pow(y[1],-1/c*log(10.0))) ;
// = A1*y[0] + A2 + A3*1/(1+ A4*pow(y[1],A5))
double A1 = (-1/t_ox);
double A2 = (a/t_ox);
double A3 = (b/t_ox);
double A4 = exp(k1/c);
double A5 = (-1/c*log(10.0));
//dydx[2] = (-1/t_ser)*y[2] + (alpha/t_ser) + (beta/t_ser)* 1/(1 + exp(k2/S)*pow(y[3],-1/S*log(10.0)) ;
// = B1*y[2] + B2 + B3* 1/(1+B4*pow(y[3],B5))
double B1 = (-1/t_ser);
double B2 = (alpha/t_ser);
double B3 = (beta/t_ser);
double B4 = exp(k2/S);
double B5 = (-1/S*log(10.0));
//dfdy[0][1] = (b*exp(k1/c)/c*log(10.0)*t_ox)*pow(y[1],-1/c*log(10.0) - 1.0)/pow(1 + exp(k1/c)*pow(y[1], 1/c*log(10.0)),2);
// = C1*(pow(y[1],C2)/pow(1 + A4*pow(y[1],A5),2);
double C1 = (b*A4/c*log(10.0)*t_ox);
double C2 = A5-1;
//double C3 = 1 + A4*pow(y[1],A5);
//dfdy[2][3] = (beta*exp(k2/S)/S*log(10.0)*t_ser)*pow(y[3],-1/Slog(10.0) - 1)/pow(1+exp(k2/S)*pow(y[3],S*log(10.0)),2)
// = D1*(pow(y[3],D2)/pow(1 + B4*pow(y[3],B5),2);
double D1 = (beta*B4/S*log(10.0)*t_ser);
double D2 = B5-1;
//double D3 = 1 + B4*pow(y[3],B5);
//==========================================================================================================================
struct rhs
{
void operator() (const Doub x, VecDoub_I &y, VecDoub_O &dydx)
{
dydx[0] = A1*y[0] + A2 + A3*1/(1 + A4*pow(y[1],A5)) ;
dydx[1] = -eta*y[1] + gamma*y[2];
dydx[2] = B1*y[2] + B2 + B3*1/(1 + B4*pow(y[3],B5)) ;
dydx[3] = phi*y[0] - Vmax*y[3]/(Km + y[3]);
}
void jacobian(const Doub x, VecDoub_I &y, VecDoub_O &dfdx,
MatDoub_O &dfdy)
{
Int n=y.size();
for (Int i=0;i<n;i++)
dfdx[i]=0.0;
dfdy[0][0] = A1;
dfdy[0][1] = C1*(pow(y[1],C2)/pow(1 + A4*pow(y[1],A5),2));
dfdy[0][2] = 0.0;
dfdy[0][3] = 0.0;
dfdy[1][0] = 0.0;
dfdy[1][1] = -eta;
dfdy[1][2] = gamma;
dfdy[1][3] = 0.0;
dfdy[2][0] = 0.0;
dfdy[2][1] = 0.0;
dfdy[2][2] = B1;
dfdy[2][3] = D1*(pow(y[3],D2)/pow(1 + B4*pow(y[3],B5),2));
dfdy[3][0] = phi;
dfdy[3][1] = 0.0;
dfdy[3][2] = 0.0;
dfdy[3][3] = -Vmax*Km/pow(Km+y[3],2);
}
};
int main()
{
ofstream outFile("data.txt"); // Opens and instantiate a data output file
const Int nvar = 4; //Number of equations to compute
const Doub atol=1.0e-5,rtol=atol, h1=2.9e-4, hmin=0.0, x1=0.0, x2=250; //temporal integration limits (x1->t1, x2->t2)
VecDoub ystart(nvar);
//Initial conditions
ystart[0]=0.5; // f_DRN(t=0) (Hz)
ystart[1]=2.8; // [Ox](t=0) (nM)
ystart[2]=5.0; // f_LHA(t=0) (Hz)
ystart[3]=1.6; // [Ser](t=0) (nM)
Output out(500); //the value inside out is the user defined number of points required!
rhs d;
Odeint<StepperRoss<rhs> > ode(ystart,x1,x2,atol,rtol,h1,hmin,out,d);
ode.integrate();
for (Int i=0;i<out.count;i++)
{
//output to console window
cout << out.xsave[i] << " " << out.ysave[0][i] << " " << out.ysave[1][i] << " " <<
out.ysave[2][i] << " " << out.ysave[3][i] << endl;
//write data to data.txt
outFile << out.xsave[i] << " " << out.ysave[0][i] << " " << out.ysave[1][i] << " " <<
out.ysave[2][i] << " " << out.ysave[3][i] << endl;
}
cout << endl << out.count << endl; //displays the number of points output by the program
outFile.close(); // Closes the data output file!
system("Pause");
}
#include <iostream>
#include "nr3.h"
#include "odeint.h"
#include "ludcmp.h"
#include "stepper.h"
#include "stepperross.h"
#include <fstream>
#include <cmath>
using namespace std;
double a = 0.3646; // parameters for equation (1) in paper
double b = 8.6971;
double k1 = 2.0732;
double c = 0.4467;
double alpha = 10.0; // parameters for equation (2) in paper
double beta = 10.0;
double k2 = 0.2041;
double S = 0.1;
double gamma = 0.77; // "rise factor" (nM)
double eta = 0.91; // "decay rate" (/s)
double t_ox = 60.0; // "time constant for orexin" (s)
double t_ser = 10.0; // "time constant for serotonin/5_HT" (s)
double phi = 33.57; //"Release per stimulus rate" (nM)
double Km = 170.0; // (nM) Km Vmax Michaelis-Menton parameters
double Vmax = 1300; // (nM/s)
//NB Input/output (3) and (4) were substiuted directly subsitituted into the equations (1) and (2) in the paper
//dydx[0] = (-1/t_ox)*y[0] + (a/t_ox) + (b/t_ox)*(1/(1 + exp(k1/c)*pow(y[1],-1/c*log(10.0))) ;
// = A1*y[0] + A2 + A3*1/(1+ A4*pow(y[1],A5))
double A1 = (-1/t_ox);
double A2 = (a/t_ox);
double A3 = (b/t_ox);
double A4 = exp(k1/c);
double A5 = (-1/c*log(10.0));
//dydx[2] = (-1/t_ser)*y[2] + (alpha/t_ser) + (beta/t_ser)* 1/(1 + exp(k2/S)*pow(y[3],-1/S*log(10.0)) ;
// = B1*y[2] + B2 + B3* 1/(1+B4*pow(y[3],B5))
double B1 = (-1/t_ser);
double B2 = (alpha/t_ser);
double B3 = (beta/t_ser);
double B4 = exp(k2/S);
double B5 = (-1/S*log(10.0));
//dfdy[0][1] = (b*exp(k1/c)/c*log(10.0)*t_ox)*pow(y[1],-1/c*log(10.0) - 1.0)/pow(1 + exp(k1/c)*pow(y[1], 1/c*log(10.0)),2);
// = C1*(pow(y[1],C2)/pow(1 + A4*pow(y[1],A5),2);
double C1 = (b*A4/c*log(10.0)*t_ox);
double C2 = A5-1;
//double C3 = 1 + A4*pow(y[1],A5);
//dfdy[2][3] = (beta*exp(k2/S)/S*log(10.0)*t_ser)*pow(y[3],-1/Slog(10.0) - 1)/pow(1+exp(k2/S)*pow(y[3],S*log(10.0)),2)
// = D1*(pow(y[3],D2)/pow(1 + B4*pow(y[3],B5),2);
double D1 = (beta*B4/S*log(10.0)*t_ser);
double D2 = B5-1;
//double D3 = 1 + B4*pow(y[3],B5);
//==========================================================================================================================
struct rhs
{
void operator() (const Doub x, VecDoub_I &y, VecDoub_O &dydx)
{
dydx[0] = A1*y[0] + A2 + A3*1/(1 + A4*pow(y[1],A5)) ;
dydx[1] = -eta*y[1] + gamma*y[2];
dydx[2] = B1*y[2] + B2 + B3*1/(1 + B4*pow(y[3],B5)) ;
dydx[3] = phi*y[0] - Vmax*y[3]/(Km + y[3]);
}
void jacobian(const Doub x, VecDoub_I &y, VecDoub_O &dfdx,
MatDoub_O &dfdy)
{
Int n=y.size();
for (Int i=0;i<n;i++)
dfdx[i]=0.0;
dfdy[0][0] = A1;
dfdy[0][1] = C1*(pow(y[1],C2)/pow(1 + A4*pow(y[1],A5),2));
dfdy[0][2] = 0.0;
dfdy[0][3] = 0.0;
dfdy[1][0] = 0.0;
dfdy[1][1] = -eta;
dfdy[1][2] = gamma;
dfdy[1][3] = 0.0;
dfdy[2][0] = 0.0;
dfdy[2][1] = 0.0;
dfdy[2][2] = B1;
dfdy[2][3] = D1*(pow(y[3],D2)/pow(1 + B4*pow(y[3],B5),2));
dfdy[3][0] = phi;
dfdy[3][1] = 0.0;
dfdy[3][2] = 0.0;
dfdy[3][3] = -Vmax*Km/pow(Km+y[3],2);
}
};
int main()
{
ofstream outFile("data.txt"); // Opens and instantiate a data output file
const Int nvar = 4; //Number of equations to compute
const Doub atol=1.0e-5,rtol=atol, h1=2.9e-4, hmin=0.0, x1=0.0, x2=250; //temporal integration limits (x1->t1, x2->t2)
VecDoub ystart(nvar);
//Initial conditions
ystart[0]=0.5; // f_DRN(t=0) (Hz)
ystart[1]=2.8; // [Ox](t=0) (nM)
ystart[2]=5.0; // f_LHA(t=0) (Hz)
ystart[3]=1.6; // [Ser](t=0) (nM)
Output out(500); //the value inside out is the user defined number of points required!
rhs d;
Odeint<StepperRoss<rhs> > ode(ystart,x1,x2,atol,rtol,h1,hmin,out,d);
ode.integrate();
for (Int i=0;i<out.count;i++)
{
//output to console window
cout << out.xsave[i] << " " << out.ysave[0][i] << " " << out.ysave[1][i] << " " <<
out.ysave[2][i] << " " << out.ysave[3][i] << endl;
//write data to data.txt
outFile << out.xsave[i] << " " << out.ysave[0][i] << " " << out.ysave[1][i] << " " <<
out.ysave[2][i] << " " << out.ysave[3][i] << endl;
}
cout << endl << out.count << endl; //displays the number of points output by the program
outFile.close(); // Closes the data output file!
system("Pause");
}
To copy to clipboard, switch view to plain text mode
Bookmarks