/* * Mathematical model predicts a critical role for osteoclast autocrine * regulation in the control of bone remodeling * * Model Status * * This CellML model runs in both OpenCell and COR to recreate * the published results (figure 2). Where parameters were not * defined in the paper they were taken from Lemaire et al. 2004. * A force of 10N was applied at day 100 and removed at day 130. * * Model Structure * * ABSTRACT: Bone is a dynamic living tissue that undergoes continuous * adaptation of its mass and structure in response to mechanical * and biological environment demands. Studies of bone adaptation * have focused on metabolic or mechanical stimulus, but mathematical * models of bone adaptation considering both, are not available * by now. In this paper, we propose a mathematical model of bone * adaptation during a remodeling cycle due to mechanical stimulus * with the introduction of osteocytes as mechanotransducers. The * model captures qualitatively very well the bone adaptation and * cell interactions during the bone remodeling. * * model diagram * * [[Image file: maldonado_2006.png]] * * Schematic representation of the basic structure of the model. * * model diagram * * [[Image file: lemaire_2004_new.png]] * * Schematic diagram of the Lemaire et al. 2004 model of which * the current model is an extension. * * The original paper reference is cited below: * * Mathematical modeling and analysis of force induced bone growth, * Maldonado S, Borchers S, Findeisen R, Allgower F, 2006, Proceedings * of the 28th IEEE EMBS Annual International Conference, 1, 3154-3157. * PubMed ID: 17947010 */ import nsrunit; // Warning: unit conversion turned off due to unit errors in 2 equation(s) unit conversion off; unit day=86400 second^1; unit pM=1E-9 meter^(-3)*mole^1; unit flux=1.1574074E-14 meter^(-3)*second^(-1)*mole^1; unit cells = fundamental; // unit picomole predefined unit picomole_cells=1E-12 mole^1*cells^1; unit picomole_per_picomole_cells=1 cells^(-1); unit picomole_per_day_picomole_cells=1.1574074E-5 second^(-1)*cells^(-1); unit per_pM=1E9 meter^3*mole^(-1); unit mm=.001 meter^1; unit mm2=1E-6 meter^2; unit mm_per_day=1.1574074E-8 meter^1*second^(-1); unit mm2_per_N=1E-6 kilogram^(-1)*meter^1*second^2; unit mm2_per_N_day=1.1574074E-11 kilogram^(-1)*meter^1*second^1; unit N_per_mm2=1E6 kilogram^1*meter^(-1)*second^(-2); unit pM_N_per_mm2=1E-3 kilogram^1*meter^(-4)*second^(-2)*mole^1; unit dyn=1E-5 kilogram^1*meter^1*second^(-2); unit mm2_day_per_dyn=8640 kilogram^(-1)*meter^1*second^3; unit first_order_rate_constant=1.1574074E-5 second^(-1); unit second_order_rate_constant=1.1574074E4 meter^3*second^(-1)*mole^(-1); unit dyn_per_mm2=10 kilogram^1*meter^(-1)*second^(-2); math main { realDomain time day; time.min=0; extern time.max; extern time.delta; real A_B(time) mm2; real r_B(time) mm; when(time=time.min) r_B=0.9912; real X_c pM; X_c=9.127e-4; real X_b pM; X_b=7.282e-4; real X_y pM; X_y=7.300e-3; real k_rB first_order_rate_constant; k_rB=1.0e0; real k_for mm_per_day; k_for=1.0e-3; real k_res mm_per_day; k_res=10.0e-3; real xc(time) pM; when(time=time.min) xc=9.127e-4; real xb(time) pM; when(time=time.min) xb=7.282e-4; real xy(time) pM; when(time=time.min) xy=7.300e-3; real X_bss pM; X_bss=7.282e-4; real X_yss pM; X_yss=7.300e-3; real k_byp first_order_rate_constant; k_byp=1.00e-1; real k_yd first_order_rate_constant; k_yd=1.00e0; real F_sti(time) pM_N_per_mm2; real k_Fs mm2_per_N; k_Fs=1.00e0; real k_y per_pM; k_y=1.00e0; real F_s(time) N_per_mm2; real F_a(time) newton; real x_no(time) pM; when(time=time.min) x_no=0.02; real X_noe flux; X_noe=0e0; real k_yno mm2_day_per_dyn; k_yno=2e4; real k_nod first_order_rate_constant; k_nod=1e3; real x_pge(time) pM; when(time=time.min) x_pge=0.01; real X_pgex flux; X_pgex=0e0; real k_pged first_order_rate_constant; k_pged=1e2; real k_nopge first_order_rate_constant; k_nopge=1e0; real k_ypge mm2_day_per_dyn; k_ypge=1e2; real x_opg(time) pM; when(time=time.min) x_opg=0.01; real k_opgd first_order_rate_constant; k_opgd=3.5e-1; real k_nopg first_order_rate_constant; k_nopg=1e1; real Io flux; Io=0.0; real K_o_p picomole_per_day_picomole_cells; K_o_p=2e5; real xr(time) pM; when(time=time.min) xr=7.734e-4; real pi_c(time) dimensionless; real x_kl(time) pM; when(time=time.min) x_kl=0.01; real pi_L(time) dimensionless; real K_l_p dimensionless; K_l_p=3e6; real K pM; K=10; real k_nokl first_order_rate_constant; k_nokl=1e2; real Il flux; Il=0; real rl flux; rl=1e3; real k1 second_order_rate_constant; k1=1e-2; real k2 first_order_rate_constant; k2=10; real k3 second_order_rate_constant; k3=5.8e-4; real k4 first_order_rate_constant; k4=1.7e-2; real ko first_order_rate_constant; ko=0.35; real pi_p dimensionless; real D_R flux; D_R=7e-4; real k_pger first_order_rate_constant; k_pger=1e-4; real D_B first_order_rate_constant; real k_B first_order_rate_constant; k_B=0.189; real D_C flux; D_C=0.189; real D_A first_order_rate_constant; D_A=0.7; real f0 dimensionless; f0=0.05; real dB first_order_rate_constant; dB=0.7; real C_s pM; C_s=5e-3; real P pM; real P_0 pM; real P_s pM; real SP flux; SP=250; real k5 second_order_rate_constant; k5=0.02; real k6 first_order_rate_constant; k6=3; real IP flux; IP=0; real kP first_order_rate_constant; kP=86; // // A_B=(r_B^2*3.141592653589793); // r_B:time=(k_for/X_b*xb+(1 mm_per_day)/X_y*xy-(k_res/X_c*xc+k_rB*r_B)); // xy:time=(k_byp*(xb-X_bss)-k_yd*(xy-X_yss)); // F_sti=(F_s*xy/(1+exp((-1)*(k_Fs*F_s+k_y*xy)))); // F_s=(F_a/A_B); F_a=(if ((time>=(100 day)) and (time<(130 day))) (1E4 dyn) else (100 dyn)); // x_no:time=(k_yno*F_sti+X_noe-k_nod*x_no); // x_pge:time=(k_ypge*F_sti+k_nopge*x_no+X_pgex-k_pged*x_pge); // x_opg:time=(K_o_p/(1 picomole_per_picomole_cells)*pi_c*xr+Io+k_nopg*x_no-k_opgd*x_opg); // x_kl:time=(rl+Il-(k_nokl*x_no+rl*((1+k1/k2*x_opg+k3/k4*K)/(K_l_p*pi_p*xb))*x_kl)); pi_L=(k3/k4*K_l_p*pi_p*xb/(1+k3*K/k4+k1/(k2*ko)*(K_o_p/(1 picomole_per_picomole_cells)/pi_p*xr+Io))*(1+Il/rl)); // xr:time=(D_R*pi_c+k_pger*x_pge-D_B/pi_c*xr); // xb:time=(D_B/pi_c*xr-k_B*xb); // xc:time=(D_C*pi_L-D_A*pi_c*xc); // D_B=(f0*dB); pi_c=((xc+f0*C_s)/(xc+C_s)); pi_p=((P+P_0)/(P+P_s)); P=(IP/kP); P_0=(SP/kP); P_s=(k6/k5); }