Skip to content
New issue

Have a question about this project? # for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “#”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? # to your account

Initial version of cell death #6

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# ABM microcarriers

This software is used for the Proof of Concept project by the [CMMC](thecmmc.org). The technical report describing the proof of concept can be found [here](https://zenodo.org/record/4264611).
This software is used for the Proof of Concept project by the [CMMC](thecmmc.org).
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It seems like you did not base this branch on the most up to date master branch. For future branches it's best to do a git pull on the master branch just before you create a new branch. For now we can just put this addition in:

Suggested change
This software is used for the Proof of Concept project by the [CMMC](thecmmc.org).
This software is used for the Proof of Concept project by the [CMMC](thecmmc.org). The technical report describing the proof of concept can be found [here](https://zenodo.org/record/4264611).


## About

Expand Down
61 changes: 40 additions & 21 deletions model_define.h
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,13 @@ typedef enum _agent_type_e {
NUM_AGENT_TYPES
} agent_type_e;

typedef enum _agent_mechanic_type_e {
MCARRIER_INERT,
CELL_A_LIVE,
CELL_A_DEATH,
NUM_MECH_TYPES,
} agent_mechanic_type_e;

typedef enum _cell_state_real_e {
CELL_MODEL_REAL_RADIUS,
CELL_MODEL_REAL_MASS,
Expand All @@ -60,9 +67,15 @@ typedef enum _cell_state_real_e {
NUM_CELL_STATE_REALS
} cell_state_real_e;

typedef enum _cell_state_int_e {
CELL_MODEL_INT_STATE,
NUM_CELL_STATE_INTS
} cell_state_int_e ;

/* ODE parameters for biomas change */
typedef enum _ode_net_GrowingCell_var_e {
ODE_NET_VAR_GROWING_CELL_BIOMASS, // dm/dt = constant
ODE_NET_VAR_STRESS_TIME, // dT/dt = heaviside
NUM_ODE_NET_VAR_GROWING_CELL
} ode_net_GrowingCell_var_e;

Expand All @@ -89,7 +102,9 @@ typedef enum _particle_extra_output_real_e {
} particle_extra_output_real_e;

typedef enum _grid_summary_real_e {
GRID_SUMMARY_REAL_MICROCARRIERS,
GRID_SUMMARY_REAL_LIVE_CELLS,
GRID_SUMMARY_REAL_DEATH,
GRID_SUMMARY_REAL_MAX_DISP,
GRID_SUMMARY_REAL_MAX_GROWRATE,
NUM_GRID_SUMMARY_REALS
Expand All @@ -109,32 +124,32 @@ typedef enum _junction_end_real_e {
} junction_end_real_e;


const REAL A_DIFFUSION_COEFF_CELLS[ NUM_AGENT_TYPES ] = { 0.0 , 0.0 };
const REAL A_AGENT_FRICIONAL_DRAG[ NUM_AGENT_TYPES ] = { 1.2e-4, 1e-5 };
const REAL A_CELL_RADIUS[ NUM_AGENT_TYPES ] = { 150.0 *0.5, 14.5 };
const REAL A_DIVISION_RADIUS[ NUM_AGENT_TYPES ] = { 0.0 , 15.0 } ;
const REAL A_MIN_CELL_RADIUS[ NUM_AGENT_TYPES ] ={ 0.0, 11.0};
const REAL A_MAX_CELL_RADIUS[ NUM_AGENT_TYPES ] ={ 0.0, 16.0};
const REAL A_MAX_CELL_VOL[NUM_AGENT_TYPES]={ 0.0, 4.0 * MY_PI * 16.0*16.0*16.0 / 3.0 };
const REAL A_DIFFUSION_COEFF_CELLS[ NUM_AGENT_TYPES ] = { 0.0, 0.0 };
const REAL A_AGENT_FRICIONAL_DRAG[ NUM_AGENT_TYPES ] = { 1.2e-4, 1e-5 };
const REAL A_CELL_RADIUS[ NUM_AGENT_TYPES ] = { 150.0 *0.5, 14.5 }; //14.5
const REAL A_DIVISION_RADIUS[ NUM_AGENT_TYPES ] = { 0.0, 15.0 };
const REAL A_MIN_CELL_RADIUS[ NUM_AGENT_TYPES ] ={ 0.0, 11.0 } ;
const REAL A_MAX_CELL_RADIUS[ NUM_AGENT_TYPES ] ={ 0.0, 16.0 };
const REAL A_MAX_CELL_VOL[NUM_AGENT_TYPES]={ 0.0, 4.0*MY_PI*16.0*16.0*16.0/3.0 };
const REAL A_MIN_CELL_VOL[NUM_AGENT_TYPES]={ 0.0, 523.6 };

// parameters of shoving adhesion
const REAL A_AGENT_SHOVING_SCALE[NUM_AGENT_TYPES] = {1.0, 1.0} ;
const REAL A_AGENT_SHOVING_LIMIT[ NUM_AGENT_TYPES ] = { 0.0 , 0.0 } ;
const REAL A_AGENT_ADHESION_S[NUM_AGENT_TYPES][NUM_AGENT_TYPES]={{0.0, 0.0}, {0.0, 0.01} };
const REAL A_AGENT_SHOVING_SCALE[ NUM_MECH_TYPES ] = {1.0, 1.0, 1.0} ;
const REAL A_AGENT_SHOVING_LIMIT[ NUM_MECH_TYPES ] = { 0.0 , 0.0, 0.0} ;
const REAL A_AGENT_ADHESION_S[NUM_MECH_TYPES][NUM_MECH_TYPES]={{0.0, 0.0, 0.0}, {0.0, 0.0,0.0},{0.0,0.0,0.0}};

// parameters of bonding forces: https://doi.org/10.1371/journal.pone.0191089
const REAL A_AGENT_BOND_S[NUM_AGENT_TYPES][NUM_AGENT_TYPES]={{0.2, 0.2}, {0.2, 0.2} };
const REAL A_AGENT_BOND_DESTROY_FACTOR[NUM_AGENT_TYPES][NUM_AGENT_TYPES] = {{0.0, 1.1}, {1.1, 1.1} };
const REAL A_AGENT_BOND_CREATE_FACTOR[NUM_AGENT_TYPES][NUM_AGENT_TYPES] = {{0.0, 1.1}, {1.1, 1.1} };
const REAL A_AGENT_STIFFNESS[NUM_AGENT_TYPES][NUM_AGENT_TYPES] = {{2.2e-3,1e-3},{1e-3,1e-3}} ;
const REAL A_AGENT_BOND_S[NUM_MECH_TYPES][NUM_MECH_TYPES]={{0.2,0.2,0.2}, {0.2,0.2,0.2},{0.2,0.2,0.2}};
const REAL A_AGENT_BOND_DESTROY_FACTOR[NUM_MECH_TYPES][NUM_MECH_TYPES] = {{0.0,1.1,1.0},{1.1,1.1,1.0},{1.0,1.0,1.0}};
const REAL A_AGENT_BOND_CREATE_FACTOR[NUM_MECH_TYPES][NUM_MECH_TYPES] = {{0.0,1.1,1.0}, {1.1,1.1,1.0},{1.0,1.0,1.0}};
const REAL A_AGENT_STIFFNESS[NUM_MECH_TYPES][NUM_MECH_TYPES] = {{2.2e-3,1e-3,1e-3},{1e-3,1e-3,1e-3},{1e-3,1e-3,1e-3}} ;


const REAL A_DENSITY_BIOMASS[ NUM_AGENT_TYPES ] = { 8.321827089772318e-17 , 8.355634512324506e-16 }; // Kg / um^3
const REAL A_DENSITY_BIOMASS[ NUM_AGENT_TYPES ] = {8.321827089772318e-17,8.355634512324506e-16 }; // Kg / um^3
const REAL DENSITY_MEDIUM = 1.0e-18 ;// Kg / um^3
const REAL A_MCARRIER_DENSITY_PER_UB = 1000 / (32.0 *32.0 * 32.0 ) ; // 2000 //0.02 ; // 0.1
const REAL INIT_CELLS_PER_MICROCARRIER = 4; //10;
const REAL A_CELL_D_MAX[ NUM_AGENT_TYPES ] = { 150.0 * 1.25, 20.0 * 1.25 };
const REAL A_MCARRIER_DENSITY_PER_UB = 1 / (32.0 *32.0 * 32.0 ) ; // 2000 //0.02 ; // 0.1
const REAL INIT_CELLS_PER_MICROCARRIER = 5; //10;
const REAL A_CELL_D_MAX[ NUM_AGENT_TYPES ] = {150.0*1.25, 20.0*1.25};


const REAL ADHESION_S = 0.01;
Expand All @@ -148,7 +163,7 @@ const REAL STEP_TIME = 1.0 ;
const REAL NUM_STATE_AND_GRID_TIME_STEPS_PER_BASELINE = 1 ;

// cell growth
const REAL DOUBLING_TIME = BASELINE_TIME_STEP_DURATION * 100000 ; // I also used 100000
const REAL DOUBLING_TIME = BASELINE_TIME_STEP_DURATION * 1000 ; // I also used 100000
const REAL ODE_CELL_GROWTH_CONSTANT = A_DENSITY_BIOMASS[1] * ( 4.0 * MY_PI * 15.0*15.0*15.0 / 3.0 ) /( 2.0 * DOUBLING_TIME ); // 4pi/3 ( R_div ^3) * density / (2 * doublingtime)
const REAL STRESS_TRESHOLD = 1e-7 ; //1e-7 ; //

Expand All @@ -170,11 +185,15 @@ const REAL DYNAMIC_VISCOSITY = 1.0e-9; // [micro N][s]/ ( [ micro m] [micro m]
const S32 AGENT_TRANSLATION_ROTATION_INTEGRATION_STEPS_PER_BASELINE_TIME_STEP = 1 ;
const REAL AGENT_TRANSLATION_ROTATION_PSEUDO_TIME_STEP_DURATION = BASELINE_TIME_STEP_DURATION / ( (REAL) AGENT_TRANSLATION_ROTATION_INTEGRATION_STEPS_PER_BASELINE_TIME_STEP );

const REAL VELOCITY_DAMPING_TEST = 1.0; //1.0e-4;
/* MODEL END */
const REAL VELOCITY_DAMPING_TEST = 0.0; // 1.0; //1.0e-4;

// gravity acceleration
const REAL STANDARD_GRAVITY = 9.80665e-6 ; // micro meter / ( s^2 )

// cell death parameters
const REAL MECH_STRESS_TRESHOLD_DEATH = 1.0 ; // 2.5e-7 ;//5e-7 ;
const REAL TIME_TO_DEATH = BASELINE_TIME_STEP_DURATION * 1000 ;

/* MODEL END */
#endif/* #ifndef __MODEL_DEFINE_H__ */

83 changes: 56 additions & 27 deletions model_routine_agent.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -73,12 +73,14 @@ void ModelRoutine::addSpAgents( const BOOL init, const VIdx& startVIdx, const VI
REAL biomass = volume_agent(A_CELL_RADIUS[AGENT_MCARRIER] ) * A_DENSITY_BIOMASS[ AGENT_MCARRIER ] ;
state.setModelReal( CELL_MODEL_REAL_MASS, biomass );
state.setModelReal( CELL_MODEL_REAL_EPS, 0.0 );
state.setModelReal( CELL_MODEL_REAL_UPTAKE_PCT, 0.0 ) ;
state.setModelReal( CELL_MODEL_REAL_UPTAKE_PCT, 1.0 ) ;
state.setModelReal( CELL_MODEL_REAL_DX, 0.0 );
state.setModelReal( CELL_MODEL_REAL_DY, 0.0 );
state.setModelReal( CELL_MODEL_REAL_DZ, 0.0 );
state.setModelReal( CELL_MODEL_REAL_STRESS, 0.0 );

state.setModelInt( CELL_MODEL_INT_STATE, MCARRIER_INERT ) ;

state.setMechIntrctBdrySphere( A_CELL_D_MAX[AGENT_MCARRIER] );

v_spAgentVIdx.push_back( vIdx );
Expand Down Expand Up @@ -142,6 +144,9 @@ void ModelRoutine::addSpAgents( const BOOL init, const VIdx& startVIdx, const VI
state_c.setModelReal( CELL_MODEL_REAL_STRESS, 0.0 );

state_c.setODEVal(0, ODE_NET_VAR_GROWING_CELL_BIOMASS, biomass );
state_c.setODEVal(0, ODE_NET_VAR_STRESS_TIME, 0.0 );

state_c.setModelInt( CELL_MODEL_INT_STATE, CELL_A_LIVE ) ;


state_c.setMechIntrctBdrySphere( A_CELL_D_MAX[ AGENT_CELL_A ] );
Expand All @@ -162,11 +167,32 @@ void ModelRoutine::addSpAgents( const BOOL init, const VIdx& startVIdx, const VI

void ModelRoutine::spAgentCRNODERHS( const S32 odeNetIdx, const VIdx& vIdx, const SpAgent& spAgent, const NbrUBEnv& nbrUBEnv, const Vector<double>& v_y, Vector<double>& v_f ) {
/* MODEL START */
// dr/dt = constant * K^2 / (K^2 + stress^2 )
REAL mech_stress = spAgent.state.getModelReal( CELL_MODEL_REAL_STRESS ) ;
REAL factor = STRESS_TRESHOLD*STRESS_TRESHOLD / ( STRESS_TRESHOLD*STRESS_TRESHOLD + mech_stress*mech_stress );
REAL r_Growth = ODE_CELL_GROWTH_CONSTANT * factor ; //* biomass;
v_f[ODE_NET_VAR_GROWING_CELL_BIOMASS]= r_Growth ;

S32 mtype = spAgent.state.getModelInt( CELL_MODEL_INT_STATE ) ;
//cout << " ODE-type " << mtype << endl;
if ( mtype == CELL_A_LIVE ) {
REAL mech_stress = spAgent.state.getModelReal( CELL_MODEL_REAL_STRESS ) ;
REAL factor = STRESS_TRESHOLD*STRESS_TRESHOLD / ( STRESS_TRESHOLD*STRESS_TRESHOLD + mech_stress*mech_stress );
REAL r_Growth = ODE_CELL_GROWTH_CONSTANT * factor ; //* biomass;
v_f[ODE_NET_VAR_GROWING_CELL_BIOMASS]= r_Growth ;

// cell death: dT/dt = Hevisaide( Stress_i - Stress_0 )
REAL heaviside = 0.0 ;
if ( mech_stress > MECH_STRESS_TRESHOLD_DEATH ) {
heaviside = 1.0;
}
v_f[ODE_NET_VAR_STRESS_TIME ]= heaviside ;
}
else if ( mtype == CELL_A_DEATH ) {
REAL r_Growth = -0.5 *ODE_CELL_GROWTH_CONSTANT ; //* biomass;
v_f[ODE_NET_VAR_GROWING_CELL_BIOMASS]= r_Growth ;
v_f[ODE_NET_VAR_STRESS_TIME ] = 0.0 ;
}
else { // error
v_f[ODE_NET_VAR_GROWING_CELL_BIOMASS]= 0.0 ;
v_f[ODE_NET_VAR_STRESS_TIME ] = 0.0 ;

}


/* MODEL END */
Expand All @@ -177,34 +203,37 @@ void ModelRoutine::spAgentCRNODERHS( const S32 odeNetIdx, const VIdx& vIdx, cons
void ModelRoutine::updateSpAgentState( const VIdx& vIdx, const JunctionData& junctionData, const VReal& vOffset, const NbrUBEnv& nbrUBEnv, SpAgentState& state/* INOUT */ ) {
/* MODEL START */

REAL uptakePct = state.getModelReal( CELL_MODEL_REAL_UPTAKE_PCT );
//REAL uptakePct = state.getModelReal( CELL_MODEL_REAL_UPTAKE_PCT );
S32 type = state.getType() ; // id of cell type



if ( ( type == AGENT_CELL_A ) && ( uptakePct > 0.0 ) ) {

//if ( ( type == AGENT_CELL_A ) && ( uptakePct > 0.0 ) ) {
if ( type == AGENT_CELL_A ) {

REAL Biomas = state.getODEVal( 0, ODE_NET_VAR_GROWING_CELL_BIOMASS );
REAL Inert = 0.0;

REAL cellVol = (Biomas + Inert)/A_DENSITY_BIOMASS[type];



if ( cellVol > A_MAX_CELL_VOL[type] ) {
cellVol = A_MAX_CELL_VOL[type] ;
Biomas = cellVol * A_DENSITY_BIOMASS[type] - Inert ;
}
CHECK( cellVol >= 0.0 );
REAL newRadius = radius_from_volume( cellVol );


state.setModelReal( CELL_MODEL_REAL_RADIUS, newRadius );
state.setModelReal( CELL_MODEL_REAL_MASS, Biomas );


//REAL ODE_NET_VAR_STRESS_TIME
REAL stress_time = state.getODEVal(0, ODE_NET_VAR_STRESS_TIME);
if ( stress_time > TIME_TO_DEATH ) {
// change the cell to cell death type
state.setModelInt( CELL_MODEL_INT_STATE, CELL_A_DEATH );
state.setModelReal( CELL_MODEL_REAL_UPTAKE_PCT, -1.0 ) ;
state.setODEVal(0, ODE_NET_VAR_STRESS_TIME, 0.0 );
}
}


/* MODEL END */
return;
Expand All @@ -229,30 +258,27 @@ void ModelRoutine::updateSpAgentBirthDeath( const VIdx& vIdx, const SpAgent& spA

S32 type = spAgent.state.getType() ;

if ( type == AGENT_CELL_A ) {
if ( type == AGENT_CELL_A ) {

//REAL rnd_num = Util::getModelRand(MODEL_RNG_UNIFORM_10PERCENT); //0.9-1.1
REAL testrad = A_DIVISION_RADIUS[type] ; //* rnd_num;
REAL cell_rad = spAgent.state.getModelReal( CELL_MODEL_REAL_RADIUS ) ;

if ((spAgent.state.getModelReal(CELL_MODEL_REAL_UPTAKE_PCT)>0.0)/* live cells */ ) {
REAL cell_rad = spAgent.state.getModelReal( CELL_MODEL_REAL_RADIUS ) ;

if( cell_rad >= testrad ) {

if ( cell_rad >= testrad ) {
if ( spAgent.state.getModelInt(CELL_MODEL_INT_STATE) == CELL_A_LIVE ) {
for( S32 i = 0 ; i < spAgent.junctionData.getNumJunctions() ; i++ ) {
const JunctionEnd& end = spAgent.junctionData.getJunctionEndRef( i );
if( end.getType() == 1 ) {
if ( end.getType() == JUNCTION_END_TYPE_MICROCARRIER ) {
divide = true;
break;
}
}

}
else if ( cell_rad <= A_MIN_CELL_RADIUS[type] ) {
disappear = true;
}
}

else if ( cell_rad <= A_MIN_CELL_RADIUS[type] ) {
//cout << "seeee " << cell_rad << " " << type << endl;
disappear = true;
}
}

/* MODEL END */
Expand Down Expand Up @@ -393,6 +419,7 @@ void ModelRoutine::divideSpAgent( const VIdx& vIdx, const JunctionData& junction
motherState.setModelReal( CELL_MODEL_REAL_MASS, mother_biomas );
motherState.setModelReal( CELL_MODEL_REAL_EPS, mother_inert );
motherState.setModelReal( CELL_MODEL_REAL_UPTAKE_PCT, 1.0 );
motherState.setModelInt( CELL_MODEL_INT_STATE, CELL_A_LIVE );

dougther_biomas = biomas - mother_biomas;
dougther_inert = inert - mother_inert;
Expand All @@ -410,10 +437,12 @@ void ModelRoutine::divideSpAgent( const VIdx& vIdx, const JunctionData& junction
daughterState.setModelReal( CELL_MODEL_REAL_MASS, dougther_biomas );
daughterState.setModelReal( CELL_MODEL_REAL_EPS, dougther_inert );
daughterState.setModelReal( CELL_MODEL_REAL_UPTAKE_PCT, 1.0 );
daughterState.setModelInt( CELL_MODEL_INT_STATE, CELL_A_LIVE );

// change the new biomass on the ODEs
motherState.setODEVal(0, ODE_NET_VAR_GROWING_CELL_BIOMASS, mother_biomas);
daughterState.setODEVal(0, ODE_NET_VAR_GROWING_CELL_BIOMASS, dougther_biomas);
daughterState.setODEVal(0, ODE_NET_VAR_STRESS_TIME, 0.0 );

// mechanical interaction range
motherState.setMechIntrctBdrySphere( A_CELL_D_MAX[ AGENT_CELL_A ] );
Expand Down
43 changes: 26 additions & 17 deletions model_routine_config.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -107,34 +107,35 @@ void ModelRoutine::updateSpAgentInfo( Vector<SpAgentInfo>& v_spAgentInfo ) {/* s

for( S32 i = 0 ; i < NUM_AGENT_TYPES; i++ ) {

ODENetInfo odeNetInfo;

/* ODE setup */
odeNetInfo.numVars = NUM_ODE_NET_VAR_GROWING_CELL;
odeNetInfo.stiff = ODE_STIFF_NORMAL;
odeNetInfo.h = 0.1;
odeNetInfo.hm = 0.01;
odeNetInfo.epsilon = 1e-6;
odeNetInfo.threshold = 1e-3;
odeNetInfo.errorThresholdVal = 0.0;
odeNetInfo.warningThresholdVal = 0.0;
odeNetInfo.setNegToZero = false;



SpAgentInfo info;

info.mechIntrctBdryType = MECH_INTRCT_BDRY_TYPE_SPHERE;
info.numStateModelReals = NUM_CELL_STATE_REALS;
info.numStateModelInts = 0;
info.numStateModelInts = NUM_CELL_STATE_INTS;
info.numStateInternalModelReals = 0;
info.numStateInternalModelInts = 0;
info.v_mechIntrctModelRealInfo.assign( NUM_CELL_MECH_REALS, modelVarInfo );
info.v_mechIntrctModelIntInfo.clear();
info.v_boolNetInfo.clear();

if ( i == AGENT_CELL_A ) {
if ( i == AGENT_CELL_A ) {

ODENetInfo odeNetInfo;

/* ODE setup cell growth */
odeNetInfo.numVars = NUM_ODE_NET_VAR_GROWING_CELL;
odeNetInfo.stiff = ODE_STIFF_NORMAL;
odeNetInfo.h = 0.1;
odeNetInfo.hm = 0.01;
odeNetInfo.epsilon = 1e-6;
odeNetInfo.threshold = 1e-3;
odeNetInfo.errorThresholdVal = 0.0;
odeNetInfo.warningThresholdVal = 0.0;
odeNetInfo.setNegToZero = true;

info.v_odeNetInfo.push_back( odeNetInfo );


} else {
info.v_odeNetInfo.clear() ;
}
Expand Down Expand Up @@ -259,10 +260,18 @@ void ModelRoutine::updateSummaryOutputInfo( Vector<SummaryOutputInfo>& v_summary
v_summaryOutputIntInfo.clear();
v_summaryOutputRealInfo.resize( NUM_GRID_SUMMARY_REALS );

info.name = "Number of microcarriers";
info.type = SUMMARY_TYPE_SUM;
v_summaryOutputRealInfo[ GRID_SUMMARY_REAL_MICROCARRIERS ] = info;

info.name = "Number of Live Cells";
info.type = SUMMARY_TYPE_SUM;
v_summaryOutputRealInfo[GRID_SUMMARY_REAL_LIVE_CELLS] = info;

info.name = "Number of Death Cells";
info.type = SUMMARY_TYPE_SUM;
v_summaryOutputRealInfo[ GRID_SUMMARY_REAL_DEATH ] = info;

info.name = "Maximum displacement";
info.type = SUMMARY_TYPE_MAX;
v_summaryOutputRealInfo[ GRID_SUMMARY_REAL_MAX_DISP ] = info;
Expand Down
Loading