Page MenuHomec4science
Files F92058537

deltaNotch4.cpp
No OneTemporary
Actions

File Metadata

Created
Sun, Nov 17, 00:26

deltaNotch4.cpp
View Options

//#include <boost/numeric/odeint.hpp>
#include "deltaNotch4.hpp"
#include <boost/math/constants/constants.hpp>
#include <cmath>
#define PI (boost::math::constants::pi<double>())
void DeltaNotch4::operator() (const state_type &y , state_type &dydt , const double t)
{
uint nbcells = lattice-> getNbCells();
dydt.assign(dydt.size(),0); // not sure if needed
for (uint c=0;c<nbcells;c++){
//vector<uint>& cell_order = lattice ->newCellOrder();
// for(vector<uint>::iterator it=cell_order.begin();it!=cell_order.end();++it){
// uint c=*it;
Cell& ce = lattice->getCell(c);
uint nbond = ce.getNbBoundaries();
double dsi = 0;
uint sind = ce.stateVectorIndex();
//cell state
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,decd;
double dj,nj,njdi;
// boundary state
double di=y[ind];
double ni=y[ind+1];
double nidi=y[ind+2];
double nidj=y[ind+3];
//extracellular state
uint ecind = b->extracellVectorIndex();
double ecd = y[ecind];
// 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;
}
// orientation factor
double ori = b->orientation();
// cout<<ori*180/PI<<" ";
// double orif = abs(sin(ori));
double orifb = 0.1+abs(cos(ori));
double orif = 0.1+abs(cos(ori));
//values from previous and next boundaries
// double np = y[b->prev()->stateVectorIndex()+1];
// double nn = y[b->next()->stateVectorIndex()+1];
double pecd = y[b->prev()->extracellVectorIndex()];
double necd = y[b->next()->extracellVectorIndex()];
//// delta
ddi= bd -gd*di; // production and degr.
ddi += ktm*njdi-ktp*di*nj*orifb; // trans binding
ddi += ks*njdi; //release after cleaving
ddi += -kcp*ni*di*orif+kcm*nidi; //// cis binding
ddi += kdep*ecd - kdem*di; // membrane binding and unbinding
//// notch
dni = bn -gn*ni; // production and degr.
dni += ktm*nidj -ktp*dj*ni*orifb; // trans binding
dni += - kcp*ni*di*orif+kcm*nidi; // cis binding
// dni += diffu*(np+nn-2*ni); //diffusion
//// cis complex
dnidi = kcp*ni*di*orif - kcm*nidi - ki*nidi;
//// trans complex
dnidj = -ktm*nidj + ktp*dj*ni*orifb - ks*nidj;
//// cleaved notch
dsi += ks*nidj;
//// extra cellular delta
//decd=0;
decd = -kdep*ecd + kdem*di-0.5*gd*ecd; // membrane binding and unbinding +degr.
decd += dem_diffu*(pecd+necd-2*ecd); // diffusion
////
dydt[ind] = ddi;
dydt[ind+1] = dni;
dydt[ind+2] = dnidi;
dydt[ind+3] = dnidj;
dydt[ecind] += decd; // this will be also updated by the next boundary
}
// cout<<endl;
dsi += -gs*si;
dydt[sind] = dsi;
}
}
uint DeltaNotch4::setStateVectorIndex(){
uint cnt=0;
visitorParamsExt vpe(&cnt, this,lattice);
lattice->forAllCellsAndBoundaries(DeltaNotch4::setCellIndex,
(void *)&vpe,
DeltaNotch4::setBoundaryIndex,
(void *)&vpe);
return cnt;
}
void DeltaNotch4::setBoundaryIndex(Boundary *bo, void *pt){
visitorParamsExt *vpe = (visitorParamsExt *) pt;
DeltaNotch2::setBoundaryIndex(bo,vpe);
uint *idx = vpe->cnt;
if(bo->extracellMatrix()==NULL){
bo->setExtracellularMatrix(vpe->lat->addExtracellularMatrix());
}
if(bo->extracellMatrix()!=NULL){
if(!bo->extracellMatrix()->indexed()){
bo->extracellMatrix()->setStateVectorIndex(*idx);
*idx += vpe->dn->extraCellStateDim();
}
}
}
void DeltaNotch4::initBoundaryState(Boundary *bo, void *pt){
DeltaNotch2::initBoundaryState(bo,pt);
state_type *y= (state_type*)pt;
uint idx = bo->extracellMatrix()->stateVectorIndex();
(*y)[idx]=0; // done redundantly
}

Event Timeline

c4science · Help