Homec4science
Diffusion tisue modeling (master) Browse

tismod/src/deltaNotch2.cpp
b7305f7b5a8bmaster

Last Change
Raw File

deltaNotch2.cpp
View Options

#include "deltaNotch2.hpp"
//using namespace std;
//using namespace boost::numeric::odeint;
//integrate<void (&)(const state_type&),state_type,double,double,double> mintegrate;
// void dummy(const state_type &x , state_type &dxdt , const double t){
// for(int i=0;i<5;i++){
// dxdt[i] = - 0.01*x[i];
// }
// }
void DeltaNotch2::operator() (const state_type &y , state_type &dydt , const double t)
{
uint nbcells = lattice-> getNbCells();
for (uint c=0;c<nbcells;c++){
Cell& ce = lattice->getCell(c);
uint nbond = ce.getNbBoundaries();
double dsi = 0;
uint sind = ce.stateVectorIndex();
double si = y[sind];
double bd = getDeltaProduction(si/nbond);
BoundaryIterator it(&ce);
for(Boundary *b=it.first();!it.isLast();b=it.next()){
uint ind = b->stateVectorIndex(); // index for cell state
Boundary *bo= b->opposite();
double ddi,dni,dnidi,dnidj;
double dj,nj,njdi;
// cell state
double ori = b->orientation();
// if(fabs(cos(ori)*sin(ori))<0.000000001){
// if((rand()/(double)RAND_MAX)<= 0.001){
// b->doT1Transition();
// }
// y[ind]=y[ind+1]=y[ind+2]=y[ind+3]=0;
// if(bo!=NULL){
// uint nind =bo->stateVectorIndex();
// y[nind]=y[nind+1]=y[nind+2]=y[nind+3]=0;
// }
// }
double di=y[ind];
double ni=y[ind+1];
double nidi=y[ind+2];
double nidj=y[ind+3];
// values from the other side of the boundary
if(bo!=NULL){
uint nind =bo->stateVectorIndex();// index for neigbouring cell
dj = y[nind];
nj = y[nind+1];
njdi = y[nind+3];
}
else{// no neigboring cell
dj=nj=njdi=0.0;
}
//values from previous and next boundaries
double np = y[b->prev()->stateVectorIndex()+1];
double nn = y[b->next()->stateVectorIndex()+1];
//// delta
ddi= bd -gd*di; // production and degr.
ddi += ktm*njdi-ktp*di*nj; // trans binding
ddi += ks*njdi; //release after cleaving
ddi += -kcp*ni*di+kcm*nidi; //// cis binding
//// notch
dni = bn -gn*ni;
dni += ktm*nidj -ktp*dj*ni; // trans binding
dni += - kcp*ni*di+kcm*nidi; // cis binding
dni += diffu*(np+nn-2*ni); //diffusion
//// cis complex
dnidi = kcp*ni*di - kcm*nidi - ki*nidi;
//// trans complex
dnidj = -ktm*nidj + ktp*dj*ni - ks*nidj;
//// cleaved notch
dsi += ks*nidj;
////
dydt[ind] = ddi;
dydt[ind+1] = dni;
dydt[ind+2] = dnidi;
dydt[ind+3] = dnidj;
}
dsi += -gs*si;
dydt[sind] = dsi;
}
}
void DeltaNotch2::writeState(){//const state_type y,const double t){
cout<<"state "<<state_dim<<": ";
for(uint i=0;i<state_dim;i++){
cout<<i<<": "<<state[i]<<"; ";
}
cout<<endl;
}
void DeltaNotch2::initState(){
void *pt = (void *)&state;
lattice->forAllCellsAndBoundaries(DeltaNotch2::initCellState,pt,DeltaNotch2::initBoundaryState,pt);
}
uint DeltaNotch2::setStateVectorIndex(){
uint cnt=0;
visitorParams vp(&cnt, this);
lattice->forAllCellsAndBoundaries(DeltaNotch2::setCellIndex,
(void *)&vp,
DeltaNotch2::setBoundaryIndex,
(void *)&vp);
return cnt;
}
void DeltaNotch2::initCellState(Cell *ce, void *pt){
state_type *y= (state_type*)pt;
// (*y)[ce->stateVectorIndex()]=0;
y->at(ce->stateVectorIndex()) = 0;
}
void DeltaNotch2::initBoundaryState(Boundary *bo, void *pt){
state_type *y= (state_type*)pt;
uint idx = bo->stateVectorIndex();
(*y)[idx]=0.9+0.2*(rand()/(double)RAND_MAX);
(*y)[idx+1]=0.9+0.2*(rand()/(double)RAND_MAX);
// (*y)[idx]=10*(rand()/(double)RAND_MAX);
//(*y)[idx+1]=10*(rand()/(double)RAND_MAX);
// (*y)[idx]=1;
//(*y)[idx+1]=1;
(*y)[idx+2]=0;
(*y)[idx+3]=0;
}
void DeltaNotch2::setCellIndex(Cell *ce, void *pt){
visitorParams *vp = (visitorParams *) pt;
uint *idx = vp->cnt;
ce->setStateVectorIndex(*idx);
*idx += vp->dn->cellStateDim();
}
void DeltaNotch2::setBoundaryIndex(Boundary *bo, void *pt){
visitorParams *vp = (visitorParams *) pt;
uint *idx = vp->cnt;
bo->setStateVectorIndex(*idx);
*idx += vp->dn->boundaryStateDim();
}
void DeltaNotch2::collectCellIdx(Cell *c,void *arg){
vector<uint> *v = (vector<uint> *)arg;
v->push_back(c->stateVectorIndex());
}
// KeepLast& DeltaNotch2::observeLast()
// {return observer;}
void KeepLast::operator()(const state_type& x, double t){
for(uint i=0;i<x.size();i++){
last[i] =x[i];
// last[i] =t;
}
}
#ifdef DELTA_NOTCH2_MAIN
int main(int argc, char* arg[]){
DeltaNotch1 dn;
state_type y = {{1.0,1.0,0.0,0.0,0.0}};
dn.setContext(1.02,0.98,0,0.3,0,0);
// size_t nsteps=integrate(dummy,y,0.0 , 10.0 , 0.1);
integrate(dn,y,0.0 , 10.0 , 0.1,writeState);
// for(uint i=0;i<5;i++){
// cout<<y[i]<<" ";
// }
// cout<<endl;
}
#endif
c4science · Help