From 2abd0f9546da45b961f5df456d9e8e90699c26b2 Mon Sep 17 00:00:00 2001 From: Boris Aguilar Date: Mon, 30 Nov 2020 18:45:18 +0000 Subject: [PATCH] Initial version of cell death --- README.md | 2 +- model_define.h | 61 ++++++++++++++++--------- model_routine_agent.cpp | 83 +++++++++++++++++++++++------------ model_routine_config.cpp | 43 +++++++++++------- model_routine_mech_intrct.cpp | 22 +++++----- model_routine_output.cpp | 22 +++++++--- run_model.xml | 2 +- 7 files changed, 153 insertions(+), 82 deletions(-) diff --git a/README.md b/README.md index b197975..b510cd5 100644 --- a/README.md +++ b/README.md @@ -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). ## About diff --git a/model_define.h b/model_define.h index f2bcec9..b9bd239 100644 --- a/model_define.h +++ b/model_define.h @@ -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, @@ -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; @@ -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 @@ -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; @@ -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 ; // @@ -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__ */ diff --git a/model_routine_agent.cpp b/model_routine_agent.cpp index 4be67ca..caac675 100644 --- a/model_routine_agent.cpp +++ b/model_routine_agent.cpp @@ -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 ); @@ -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 ] ); @@ -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& v_y, Vector& 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 */ @@ -177,20 +203,16 @@ 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] ; @@ -198,13 +220,20 @@ void ModelRoutine::updateSpAgentState( const VIdx& vIdx, const JunctionData& jun } 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; @@ -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 */ @@ -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; @@ -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 ] ); diff --git a/model_routine_config.cpp b/model_routine_config.cpp index 896d346..d7c1de1 100644 --- a/model_routine_config.cpp +++ b/model_routine_config.cpp @@ -107,34 +107,35 @@ void ModelRoutine::updateSpAgentInfo( Vector& 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() ; } @@ -259,10 +260,18 @@ void ModelRoutine::updateSummaryOutputInfo( Vector& 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; diff --git a/model_routine_mech_intrct.cpp b/model_routine_mech_intrct.cpp index cb66fd4..8004a6b 100644 --- a/model_routine_mech_intrct.cpp +++ b/model_routine_mech_intrct.cpp @@ -35,7 +35,6 @@ void ModelRoutine::initJunctionSpAgent( const VIdx& vIdx0, const SpAgent& spAgen REAL dist_threshold = factor * (R0 + R1); if ( A_AGENT_BOND_S[type0][type1] > 0.0 ){ - //if ( type0 != type1 ) { if ( dist < dist_threshold ) { link = true; @@ -68,19 +67,22 @@ void ModelRoutine::computeMechIntrctSpAgent( const S32 iter, const VIdx& vIdx0, S32 type0 = spAgent0.state.getType(); S32 type1 = spAgent1.state.getType(); - REAL R0 = A_AGENT_SHOVING_SCALE[type0]*spAgent0.state.getModelReal( CELL_MODEL_REAL_RADIUS ) ; - REAL R1 = A_AGENT_SHOVING_SCALE[type1]*spAgent1.state.getModelReal( CELL_MODEL_REAL_RADIUS ) ; + S32 mtype0 = spAgent0.state.getModelInt( CELL_MODEL_INT_STATE ) ; + S32 mtype1 = spAgent1.state.getModelInt( CELL_MODEL_INT_STATE ) ; + + REAL R0 = A_AGENT_SHOVING_SCALE[mtype0]*spAgent0.state.getModelReal( CELL_MODEL_REAL_RADIUS ) ; + REAL R1 = A_AGENT_SHOVING_SCALE[mtype1]*spAgent1.state.getModelReal( CELL_MODEL_REAL_RADIUS ) ; REAL dist_threshold = R0 + R1; - REAL D = R0 + R1 - 0.5*A_AGENT_SHOVING_LIMIT[type0] - 0.5*A_AGENT_SHOVING_LIMIT[type1]; + REAL D = R0 + R1 - 0.5*A_AGENT_SHOVING_LIMIT[mtype0] - 0.5*A_AGENT_SHOVING_LIMIT[mtype1]; REAL mag = 0.0; REAL stress = 0.0 ; REAL xij = D - dist ; - if ( A_AGENT_BOND_S[type0][type1] > 0.0 ){ - REAL sij = A_AGENT_BOND_S[type0][type1] ; + if ( A_AGENT_BOND_S[mtype0][mtype1] > 0.0 ){ + REAL sij = A_AGENT_BOND_S[mtype0][mtype1] ; if(spAgent0.junctionData.isLinked(spAgent1.junctionData) == true) { - if( dist > A_AGENT_BOND_DESTROY_FACTOR[type0][type1]* dist_threshold ) { + if( dist > A_AGENT_BOND_DESTROY_FACTOR[mtype0][mtype1]* dist_threshold ) { unlink = true;/* break junction */ //cout << " broken juntion " << endl; } @@ -105,7 +107,7 @@ void ModelRoutine::computeMechIntrctSpAgent( const S32 iter, const VIdx& vIdx0, } } else {/* no junction */ - if( dist < A_AGENT_BOND_CREATE_FACTOR[type0][type1]*dist_threshold ) { + if( dist < A_AGENT_BOND_CREATE_FACTOR[mtype0][mtype1]*dist_threshold ) { link = true;/* form junction */ @@ -131,8 +133,8 @@ void ModelRoutine::computeMechIntrctSpAgent( const S32 iter, const VIdx& vIdx0, //if ( type0 == type1 ) { // cout << "force " << mag << " xij " << xij << endl; //} - mag = A_AGENT_STIFFNESS[type0][type1] * mag ; - stress = A_AGENT_STIFFNESS[type0][type1] * stress ; + mag = A_AGENT_STIFFNESS[mtype0][mtype1] * mag ; + stress = A_AGENT_STIFFNESS[mtype0][mtype1] * stress ; mechIntrctData0.setModelReal(CELL_MECH_REAL_FORCE_X,vDir[0]*mag); diff --git a/model_routine_output.cpp b/model_routine_output.cpp index f405947..008b762 100644 --- a/model_routine_output.cpp +++ b/model_routine_output.cpp @@ -24,7 +24,7 @@ using namespace std; void ModelRoutine::updateSpAgentOutput( const VIdx& vIdx, const SpAgent& spAgent, REAL& color, Vector& v_extraReal, Vector& v_extraVReal ) { /* MODEL START */ - color = spAgent.state.getType(); + color = spAgent.state.getModelInt( CELL_MODEL_INT_STATE ); CHECK( v_extraReal.size() == NUM_PARTICLE_EXTRA_OUTPUT_REALS ); REAL radius = spAgent.state.getModelReal( CELL_MODEL_REAL_RADIUS ); v_extraReal[PARTICLE_EXTRA_OUTPUT_REAL_RADIUS] = radius; @@ -66,15 +66,25 @@ void ModelRoutine::updateSummaryVar( const VIdx& vIdx, const NbrUBAgentData& nbr const UBAgentData& ubAgentData = *( nbrUBAgentData.getConstPtr( 0, 0, 0 ) ); - REAL count = 0.0; + REAL count_micro = 0.0; + REAL count_live = 0.0; + REAL count_death = 0.0; REAL max_disp = 0.0; REAL max_fact = -1.0 ; /* Count the number of cells placed in the Simulation Domain */ for (S32 i = 0 ; i < ( S32 )ubAgentData.v_spAgent.size() ; i++ ) { + + S32 mtype = ubAgentData.v_spAgent[i].state.getModelInt( CELL_MODEL_INT_STATE ); + + if ( mtype == MCARRIER_INERT ) + count_micro += 1.0 ; + else if ( mtype == CELL_A_LIVE ) + count_live += 1.0; + else if ( mtype == CELL_A_DEATH ) + count_death += 1.0; - count += 1.0 ; REAL dx = FABS( ubAgentData.v_spAgent[i].state.getModelReal( CELL_MODEL_REAL_DX )) ; if ( dx > max_disp ) max_disp = dx ; @@ -90,9 +100,11 @@ void ModelRoutine::updateSummaryVar( const VIdx& vIdx, const NbrUBAgentData& nbr } - + /* GRID_SUMMARY_REAL_LIVE_CELLS is set in model_routine_config.cpp */ - v_realVal[GRID_SUMMARY_REAL_LIVE_CELLS] = count; + v_realVal[ GRID_SUMMARY_REAL_MICROCARRIERS ] = count_micro; + v_realVal[ GRID_SUMMARY_REAL_LIVE_CELLS ] = count_live; + v_realVal[ GRID_SUMMARY_REAL_DEATH ] = count_death; v_realVal[GRID_SUMMARY_REAL_MAX_DISP ] = max_disp * BASELINE_TIME_STEP_DURATION ; v_realVal[GRID_SUMMARY_REAL_MAX_GROWRATE ] = max_fact ; /* MODEL END */ diff --git a/run_model.xml b/run_model.xml index 00741d4..75d246e 100644 --- a/run_model.xml +++ b/run_model.xml @@ -8,7 +8,7 @@ start: starting baseline time step, required, integer, [0,) end: ending baseline time step, required, integer, [start + 1,) --> - +