Repository navigation
Expand file tree
/
Copy pathcml_mrgsolve_model.R
More file actions
714 lines (621 loc) · 28.3 KB
/
Copy pathcml_mrgsolve_model.R
File metadata and controls
714 lines (621 loc) · 28.3 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
##############################################################################
# CML QSP Model — mrgsolve ODE Implementation
# Chronic Myeloid Leukemia: BCR-ABL Kinetics, TKI PK/PD, Resistance, TFR
#
# Framework: Michor et al. (2005) Nature 435:1267-1270 (leukemic hierarchy)
# + TKI PK/PD (IRIS, ENESTnd, DASISION, ASCEMBL calibration)
# + Sensitive/Resistant clone dynamics (T315I)
#
# Compartments (22 ODEs):
# Normal hematopoiesis: y0 (HSC), y1 (CMP), y2 (GMP), y3 (Neutrophil)
# CML leukemic: x0q (LSC quiescent), x0p (LSC cycling),
# x1 (LPC), x2 (LPC diff), x3 (WBC/blast)
# Sensitive clone: xs (TKI-sensitive LSC), xsp (sensitive progenitors)
# Resistant clone: xr (T315I LSC), xrp (resistant progenitors)
# Imatinib PK: Gut_imt, Cp_imt, Cic_imt (plasma + intracellular)
# Dasatinib PK: Cp_das, Cic_das
# Nilotinib PK: Gut_nil, Cp_nil, Cic_nil
# Asciminib PK: Gut_asc, Cp_asc
# BCR-ABL IS%: computed from x3/(x3+y3) ratio
#
# Author: Generated by Claude Code (CCR) — 2026-06-20
# References: Michor 2005; Roeder 2006; IRIS 2003; ENESTnd 2010;
# DASISION 2010; ASCEMBL 2021; OPTIC 2022
##############################################################################
library(mrgsolve)
library(dplyr)
library(ggplot2)
library(tidyr)
# ============================================================
# MODEL CODE
# ============================================================
cml_code <- '
$PROB
CML QSP Model — BCR-ABL Hierarchy + TKI PK/PD + Resistance
Calibrated to IRIS, ENESTnd, DASISION, ASCEMBL clinical trials.
$PARAM @annotated
// --- Normal Hematopoiesis Parameters ---
ay0 : 0.008 : HSC self-renewal rate (Michor 2005 calibrated)
by0 : 0.004 : HSC differentiation rate
dy0 : 0.003 : HSC death rate
ay1 : 0.0 : CMP self-renewal (transit amplifying)
by1 : 0.05 : CMP differentiation rate
dy1 : 0.02 : CMP death rate
by2 : 0.2 : GMP differentiation rate
dy2 : 0.05 : GMP death rate
dy3 : 0.10 : Mature cell (neutrophil) death rate
// Amplification factors (normal)
c1 : 10.0 : HSC→CMP amplification
c2 : 8.0 : CMP→GMP amplification
c3 : 16.0 : GMP→Neutrophil amplification
// --- CML Leukemic Hierarchy Parameters ---
// LSC (x0q = quiescent, x0p = cycling)
a0 : 0.04 : LSC self-renewal rate (BCR-ABL driven, > normal)
b0 : 0.01 : LSC differentiation rate
d0 : 0.003 : LSC death rate (baseline, low due to BCR-ABL)
q_act : 0.005 : Quiescent LSC activation rate (G0→cycling)
q_dorm: 0.02 : LSC quiescence induction rate
// LPC proliferation
a1 : 0.03 : LPC1 self-renewal
b1 : 0.08 : LPC1 differentiation
d1 : 0.015 : LPC1 death rate
b2 : 0.15 : LPC2 differentiation
d2 : 0.04 : LPC2 death rate
d3 : 0.12 : Leukemic blast death rate
// Amplification (leukemic, higher than normal)
lc1 : 15.0 : LSC→LPC1 amplification
lc2 : 10.0 : LPC1→LPC2 amplification
lc3 : 20.0 : LPC2→WBC amplification
// --- Sensitive Clone Parameters ---
as0 : 0.035 : Sensitive LSC self-renewal
bs0 : 0.012 : Sensitive LSC differentiation
ds0 : 0.004 : Sensitive LSC death
asp : 0.04 : Sensitive progenitor proliferation
dsp : 0.02 : Sensitive progenitor death
// --- Resistant Clone (T315I) Parameters ---
ar0 : 0.032 : Resistant LSC self-renewal (slight fitness cost)
br0 : 0.010 : Resistant LSC differentiation
dr0 : 0.005 : Resistant LSC death
arp : 0.038 : Resistant progenitor proliferation
drp : 0.022 : Resistant progenitor death
mu : 1e-8 : Mutation rate per cell division (WT→T315I)
// --- Imatinib PK Parameters (IRIS trial) ---
// 400 mg/day standard; 800 mg/day high-dose
ka_imt : 0.8 : Absorption rate constant (h^-1), Tmax~2h
F_imt : 0.98 : Oral bioavailability
Vd_imt : 435.0 : Volume of distribution (L), Cp,ss ~1-2 µM
CL_imt : 14.0 : Total clearance (L/h), t1/2~18h
Kin_imt : 0.5 : OCT1-mediated influx rate constant (h^-1)
Kout_imt: 0.3 : Intracellular efflux (ABCB1) rate constant (h^-1)
// Dose-response: IC50=0.1-0.6 µM for BCR-ABL p210
IC50_imt : 0.25 : IC50 imatinib (µM, intracellular BCR-ABL inhibition)
Emax_imt : 0.90 : Maximum effect (90% kill of sensitive cells)
Hill_imt : 1.2 : Hill coefficient
// --- Dasatinib PK Parameters (DASISION trial) ---
ka_das : 3.5 : Absorption rate (h^-1), Tmax~0.5h
F_das : 0.24 : Bioavailability (14-34%)
Vd_das : 2505.0 : Volume of distribution (L), highly distributed
CL_das : 220.0 : Total clearance (L/h), t1/2~5-6h
// Intracellular: passive, no OCT1 dependence
Kic_das : 0.8 : Intracellular accumulation rate
Kout_das: 0.4 : Intracellular efflux
IC50_das : 0.025 : IC50 dasatinib (µM, ~10x more potent than imatinib)
Emax_das : 0.95 : Maximum effect
Hill_das : 1.1 : Hill coefficient
// --- Nilotinib PK Parameters (ENESTnd trial) ---
// 300 mg BID; fasting required
ka_nil : 0.6 : Absorption rate (h^-1), slower (fasting)
F_nil : 0.35 : Bioavailability (fasting; food ↑4-fold but CL ↑)
Vd_nil : 3100.0 : Volume of distribution (L)
CL_nil : 42.0 : Total clearance (L/h), t1/2~17h
Kic_nil : 0.6 : Intracellular accumulation
Kout_nil: 0.35 : Intracellular efflux
IC50_nil : 0.02 : IC50 nilotinib (µM, ~20x more potent than imatinib)
Emax_nil : 0.93 : Maximum effect
Hill_nil : 1.15 : Hill coefficient
// --- Asciminib PK Parameters (ASCEMBL, OPTIC trials) ---
// 40 mg BID standard; 200 mg BID for T315I
ka_asc : 1.2 : Absorption rate (h^-1)
F_asc : 0.72 : Bioavailability (~72%)
Vd_asc : 95.0 : Volume of distribution (L)
CL_asc : 7.5 : Total clearance (L/h), t1/2~8.7h
IC50_asc : 0.0005 : IC50 asciminib (µM, 0.5 nM, very potent STAMP)
IC50_asc_r : 0.004 : IC50 asciminib vs T315I (200mg BID needed, µM)
Emax_asc : 0.90 : Maximum effect (sensitive)
Emax_asc_r : 0.82 : Maximum effect vs T315I
Hill_asc : 1.3 : Hill coefficient
// --- Tumor suppression / competition ---
Kcomp : 1e-9 : Competition coefficient (density dependence)
K_total : 1e12 : Carrying capacity (total cells in marrow)
// --- Drug switches (0=off, 1=on) ---
dose_imt : 0 : Imatinib dose (mg/day)
dose_das : 0 : Dasatinib dose (mg/day)
dose_nil : 0 : Nilotinib dose (mg/day, BID; specify daily total)
dose_asc : 0 : Asciminib dose (mg/day, BID; specify daily total)
// --- Immune parameters ---
kNK_kill : 0.001 : NK cell kill rate of leukemic cells
kCTL_kill: 0.0005 : CTL kill rate of leukemic cells
NK_base : 1e8 : Baseline NK cell count
CTL_base : 5e7 : Baseline CTL count
NK_imt_effect : 0.3 : Imatinib immune restoration (% NK increase)
$CMT @annotated
// Normal hematopoiesis
y0 : HSC (normal hematopoietic stem cells)
y1 : CMP (common myeloid progenitor, normal)
y2 : GMP (granulocyte-macrophage progenitor, normal)
y3 : Neutrophils_normal (mature blood cells)
// Leukemic hierarchy (total CML)
x0q : LSC_quiescent (BCR-ABL+ quiescent, TKI-refractory)
x0p : LSC_cycling (BCR-ABL+ proliferating, TKI-sensitive)
x1 : LPC1 (leukemic progenitor committed)
x2 : LPC2 (leukemic progenitor differentiated)
x3 : WBC_leukemic (mature CML cells / blasts)
// Sensitive clone
xs : LSC_sensitive (TKI-sensitive, BCR-ABL WT)
xsp : LPC_sensitive (sensitive progenitors)
// Resistant clone (T315I-like)
xr : LSC_resistant (T315I/resistant LSC)
xrp : LPC_resistant (resistant progenitors)
// Imatinib PK
Gut_imt : Imatinib_gut
Cp_imt : Imatinib_plasma_uM
Cic_imt : Imatinib_intracellular_uM
// Dasatinib PK
Cp_das : Dasatinib_plasma_uM
Cic_das : Dasatinib_intracellular_uM
// Nilotinib PK
Gut_nil : Nilotinib_gut
Cp_nil : Nilotinib_plasma_uM
Cic_nil : Nilotinib_intracellular_uM
// Asciminib PK
Gut_asc : Asciminib_gut
Cp_asc : Asciminib_plasma_uM
$INIT @annotated
// Steady-state normal hematopoiesis (Michor 2005 calibrated)
y0 = 1e4 : HSC initial count
y1 = 5e4 : CMP initial count
y2 = 4e5 : GMP initial count
y3 = 4e9 : Normal WBC (4×10^9 / L × ~1L BM equivalent)
// CML at diagnosis (CP phase): ~10^12 leukemic cells total
x0q = 1e4 : Quiescent LSC at diagnosis
x0p = 2e3 : Cycling LSC at diagnosis
x1 = 1e7 : LPC1 at diagnosis
x2 = 5e7 : LPC2 at diagnosis
x3 = 5e10 : Leukemic WBC at diagnosis (WBC ~80×10^9/L)
xs = 8e3 : Sensitive LSC (majority)
xsp = 8e6 : Sensitive LPC
xr = 10 : Resistant LSC (rare at diagnosis, ~10 cells)
xrp = 100 : Resistant LPC
// PK: all zero at t=0
Gut_imt = 0
Cp_imt = 0
Cic_imt = 0
Cp_das = 0
Cic_das = 0
Gut_nil = 0
Cp_nil = 0
Cic_nil = 0
Gut_asc = 0
Cp_asc = 0
$ODE
// ============================================================
// Helper: MW conversions for dose (mg) → µmol
// Imatinib MW = 493.6 g/mol; Dasatinib = 488.0; Nilotinib = 529.5; Asciminib = 449.5
// Dose entered as mg/day; rate into gut = Dose/24 mg/h (continuous infusion approx)
// ============================================================
double dose_rate_imt = dose_imt / 24.0 / 493.6 * 1000.0; // µmol/h
double dose_rate_das = dose_das / 24.0 / 488.0 * 1000.0;
double dose_rate_nil = dose_nil / 24.0 / 529.5 * 1000.0;
double dose_rate_asc = dose_asc / 24.0 / 449.5 * 1000.0;
// ============================================================
// PK ODEs
// ============================================================
// --- Imatinib PK ---
dxdt_Gut_imt = F_imt * dose_rate_imt - ka_imt * Gut_imt;
dxdt_Cp_imt = ka_imt * Gut_imt / Vd_imt - (CL_imt / Vd_imt) * Cp_imt;
dxdt_Cic_imt = Kin_imt * Cp_imt - Kout_imt * Cic_imt;
// --- Dasatinib PK ---
// Continuous approximation (rapid Tmax, large Vd)
dxdt_Cp_das = F_das * dose_rate_das / Vd_das - (CL_das / Vd_das) * Cp_das;
dxdt_Cic_das = Kic_das * Cp_das - Kout_das * Cic_das;
// --- Nilotinib PK ---
dxdt_Gut_nil = F_nil * dose_rate_nil - ka_nil * Gut_nil;
dxdt_Cp_nil = ka_nil * Gut_nil / Vd_nil - (CL_nil / Vd_nil) * Cp_nil;
dxdt_Cic_nil = Kic_nil * Cp_nil - Kout_nil * Cic_nil;
// --- Asciminib PK ---
dxdt_Gut_asc = F_asc * dose_rate_asc - ka_asc * Gut_asc;
dxdt_Cp_asc = ka_asc * Gut_asc / Vd_asc - (CL_asc / Vd_asc) * Cp_asc;
// ============================================================
// PD: Drug Effect Functions (Emax model, Hill equation)
// Effect = fraction of BCR-ABL-driven proliferation inhibited
// ============================================================
// Imatinib effect on sensitive (WT BCR-ABL) cells
double Cimt_eff = (Cic_imt > 0) ? Cic_imt : 0;
double E_imt = Emax_imt * pow(Cimt_eff, Hill_imt) /
(pow(IC50_imt, Hill_imt) + pow(Cimt_eff, Hill_imt));
// Dasatinib effect on sensitive cells
double Cdas_eff = (Cic_das > 0) ? Cic_das : 0;
double E_das = Emax_das * pow(Cdas_eff, Hill_das) /
(pow(IC50_das, Hill_das) + pow(Cdas_eff, Hill_das));
// Nilotinib effect on sensitive cells
double Cnil_eff = (Cic_nil > 0) ? Cic_nil : 0;
double E_nil = Emax_nil * pow(Cnil_eff, Hill_nil) /
(pow(IC50_nil, Hill_nil) + pow(Cnil_eff, Hill_nil));
// Asciminib effect on sensitive cells (40 mg BID)
double Casc_eff = (Cp_asc > 0) ? Cp_asc : 0;
double E_asc = Emax_asc * pow(Casc_eff, Hill_asc) /
(pow(IC50_asc, Hill_asc) + pow(Casc_eff, Hill_asc));
// Asciminib effect on T315I resistant cells (200 mg BID — higher IC50)
double E_asc_r = Emax_asc_r * pow(Casc_eff, Hill_asc) /
(pow(IC50_asc_r, Hill_asc) + pow(Casc_eff, Hill_asc));
// Combined TKI effect (additive, capped at 0.98)
double E_total_sens = fmin(0.98, E_imt + E_das + E_nil + E_asc);
double E_total_res = fmin(0.98, E_asc_r); // Only asciminib (high dose) hits T315I
// ============================================================
// Density-dependent competition (Lotka-Volterra type)
// ============================================================
double Ntotal = x0q + x0p + x1 + x2 + x3 + xs + xsp + xr + xrp +
y0 + y1 + y2 + y3;
double suppress = fmax(0.0, 1.0 - Kcomp * Ntotal);
// ============================================================
// Normal Hematopoiesis ODEs
// ============================================================
// y0: HSC — self-renewal + input - differentiation - death
dxdt_y0 = ay0 * y0 + suppress * 0.0 - by0 * y0 - dy0 * y0;
dxdt_y1 = c1 * by0 * y0 - by1 * y1 - dy1 * y1;
dxdt_y2 = c2 * by1 * y1 - by2 * y2 - dy2 * y2;
dxdt_y3 = c3 * by2 * y2 - dy3 * y3;
// ============================================================
// Leukemic Hierarchy ODEs (total CML pool: xs+xr tracking separately)
// ============================================================
// LSC quiescent: activated by TKI pressure → quiescence; killed by NK
double q_in = q_dorm * x0p; // cycling→quiescent (increases under TKI)
double q_out = q_act * x0q; // quiescent→cycling
double NK_kill_x0q = kNK_kill * NK_base * x0q / (x0q + 1e3);
dxdt_x0q = q_in - q_out - d0 * x0q - NK_kill_x0q;
// LSC cycling: self-renewal + activation - quiescence - differentiation - death - TKI kill - NK
double TKI_kill_x0p = E_total_sens * x0p;
double NK_kill_x0p = kNK_kill * NK_base * x0p / (x0p + 1e3);
dxdt_x0p = a0 * x0p + q_out - q_in - b0 * x0p - d0 * x0p
- TKI_kill_x0p - NK_kill_x0p + mu * a0 * xs; // mutation from xs
// LPC1
double TKI_kill_x1 = E_total_sens * x1 * 0.6; // progenitors less sensitive
dxdt_x1 = lc1 * b0 * x0p - b1 * x1 - d1 * x1 - TKI_kill_x1;
// LPC2
double TKI_kill_x2 = E_total_sens * x2 * 0.8;
dxdt_x2 = lc2 * b1 * x1 - b2 * x2 - d2 * x2 - TKI_kill_x2;
// Leukemic WBC/blasts
double TKI_kill_x3 = E_total_sens * x3 * 1.0; // most sensitive (post-mitotic)
double NK_kill_x3 = kNK_kill * NK_base * x3 / (x3 + 1e8);
double CTL_kill_x3 = kCTL_kill * CTL_base * x3 / (x3 + 1e8);
dxdt_x3 = lc3 * b2 * x2 - d3 * x3 - TKI_kill_x3 - NK_kill_x3 - CTL_kill_x3;
// ============================================================
// Sensitive Clone ODEs
// ============================================================
double TKI_kill_xs = E_total_sens * xs;
dxdt_xs = as0 * xs - bs0 * xs - ds0 * xs - TKI_kill_xs - mu * as0 * xs;
dxdt_xsp = lc1 * bs0 * xs + asp * xsp - dsp * xsp - E_total_sens * xsp * 0.7;
// ============================================================
// Resistant Clone (T315I) ODEs
// ============================================================
double TKI_kill_xr = E_total_res * xr; // only asciminib 200mg hits this
double TKI_kill_xrp = E_total_res * xrp * 0.7;
dxdt_xr = ar0 * xr - br0 * xr - dr0 * xr - TKI_kill_xr;
dxdt_xrp = lc1 * br0 * xr + arp * xrp - drp * xrp - TKI_kill_xrp;
$TABLE
// BCR-ABL IS% calculation
// IS% = CML cells / (CML + normal cells) × 100
// Approximate: total CML WBC / (total CML WBC + normal neutrophils)
double CML_total = x3 + xrp + xsp;
double Norm_total = y3;
double BCRPCT = 100.0 * CML_total / fmax(CML_total + Norm_total, 1.0);
// Log10 BCR-ABL IS
double LOG_IS = (BCRPCT > 0.0001) ? log10(BCRPCT) : -5.0;
// Molecular response milestones
double CHR_reached = (BCRPCT < 10.0) ? 1.0 : 0.0;
double MMR_reached = (BCRPCT < 0.1) ? 1.0 : 0.0;
double MR4_reached = (BCRPCT < 0.01) ? 1.0 : 0.0;
double MR45_reached = (BCRPCT < 0.0032) ? 1.0 : 0.0;
// WBC count (10^9/L, assuming 5L blood volume)
double WBC_total_blood = (x3 + Norm_total) / 5e9; // cells/L / (cells/count unit)
// PK summary
double Cp_active = Cp_imt + Cp_das + Cp_nil + Cp_asc;
// Resistant fraction
double Resist_frac = (xr + xrp) / fmax(CML_total, 1.0);
capture BCRPCT = BCRPCT;
capture LOG_IS = LOG_IS;
capture CHR = CHR_reached;
capture MMR = MMR_reached;
capture MR4 = MR4_reached;
capture MR45 = MR45_reached;
capture WBC = WBC_total_blood;
capture LSC_total = x0q + x0p + xs + xr;
capture LSC_quies = x0q / fmax(x0q + x0p, 1.0); // quiescent fraction
capture Resist_frac = Resist_frac;
capture E_imt_out = E_imt;
capture E_das_out = E_das;
capture E_nil_out = E_nil;
capture E_asc_out = E_asc;
$CAPTURE BCRPCT LOG_IS CHR MMR MR4 MR45 WBC LSC_total LSC_quies
Resist_frac E_imt_out E_das_out E_nil_out E_asc_out
Cp_imt Cic_imt Cp_das Cic_das Cp_nil Cic_nil Cp_asc
'
# ============================================================
# Compile model
# ============================================================
mod <- mcode("cml_qsp", cml_code)
cat("Model compiled successfully. Compartments:", length(mod@cmts), "\n")
cat("Parameters:", length(param(mod)), "\n")
# ============================================================
# SIMULATION SCENARIOS
# ============================================================
# Common simulation time: 10 years = 87,600 hours
t_end_yr <- 10 # years
t_end_hr <- t_end_yr * 365.25 * 24 # hours
dt <- 24 # output every 24h (daily)
times <- seq(0, t_end_hr, by = dt)
# Convert to months for display
to_months <- function(t_hr) t_hr / (365.25 * 24 / 12)
#-----------------------------------------------------------------
# SCENARIO 1: Untreated CML progression (CP → potential AP/BC)
#-----------------------------------------------------------------
cat("\n--- Scenario 1: Untreated CML ---\n")
scen1 <- mod %>%
param(dose_imt = 0, dose_das = 0, dose_nil = 0, dose_asc = 0) %>%
mrgsim(end = t_end_hr, delta = dt) %>%
as.data.frame() %>%
mutate(scenario = "Untreated", time_months = to_months(time))
cat("Final BCR-ABL IS% (untreated, 10yr):", tail(scen1$BCRPCT, 1), "%\n")
#-----------------------------------------------------------------
# SCENARIO 2: Imatinib 400 mg/day (IRIS trial — standard of care 2001)
#-----------------------------------------------------------------
cat("\n--- Scenario 2: Imatinib 400mg/day (IRIS) ---\n")
scen2 <- mod %>%
param(dose_imt = 400, dose_das = 0, dose_nil = 0, dose_asc = 0) %>%
mrgsim(end = t_end_hr, delta = dt) %>%
as.data.frame() %>%
mutate(scenario = "Imatinib 400mg/day", time_months = to_months(time))
# Check IRIS-calibrated milestones:
# CHR @ 3mo: ~95%; CCyR @ 12mo: ~73%; MMR @ 18mo: ~39%
mmr_3mo <- scen2$MMR[scen2$time_months >= 3][1]
mmr_12mo <- scen2$MMR[scen2$time_months >= 12][1]
mmr_18mo <- scen2$MMR[scen2$time_months >= 18][1]
cat("MMR achieved @ 3mo:", mmr_3mo, "| 12mo:", mmr_12mo, "| 18mo:", mmr_18mo, "\n")
cat("BCR-ABL IS% @ 12mo:", round(scen2$BCRPCT[scen2$time_months >= 12][1], 3), "%\n")
#-----------------------------------------------------------------
# SCENARIO 3: Dasatinib 100 mg/day (DASISION trial)
# Faster and deeper responses vs imatinib
#-----------------------------------------------------------------
cat("\n--- Scenario 3: Dasatinib 100mg/day (DASISION) ---\n")
scen3 <- mod %>%
param(dose_imt = 0, dose_das = 100, dose_nil = 0, dose_asc = 0) %>%
mrgsim(end = t_end_hr, delta = dt) %>%
as.data.frame() %>%
mutate(scenario = "Dasatinib 100mg/day", time_months = to_months(time))
mmr_12mo_das <- scen3$MMR[scen3$time_months >= 12][1]
cat("MMR @ 12mo (Dasatinib):", mmr_12mo_das, "\n")
cat("BCR-ABL IS% @ 12mo:", round(scen3$BCRPCT[scen3$time_months >= 12][1], 3), "%\n")
#-----------------------------------------------------------------
# SCENARIO 4: Nilotinib 300mg BID (ENESTnd trial)
# Superior MMR rates vs imatinib
#-----------------------------------------------------------------
cat("\n--- Scenario 4: Nilotinib 300mg BID (ENESTnd) ---\n")
scen4 <- mod %>%
param(dose_imt = 0, dose_das = 0, dose_nil = 600, dose_asc = 0) %>%
# 300mg BID = 600mg/day total
mrgsim(end = t_end_hr, delta = dt) %>%
as.data.frame() %>%
mutate(scenario = "Nilotinib 300mg BID", time_months = to_months(time))
mmr_12mo_nil <- scen4$MMR[scen4$time_months >= 12][1]
mr4_12mo_nil <- scen4$MR4[scen4$time_months >= 12][1]
cat("MMR @ 12mo (Nilotinib):", mmr_12mo_nil, "\n")
cat("MR4 @ 12mo (Nilotinib):", mr4_12mo_nil, "\n")
cat("BCR-ABL IS% @ 12mo:", round(scen4$BCRPCT[scen4$time_months >= 12][1], 3), "%\n")
#-----------------------------------------------------------------
# SCENARIO 5: Asciminib 40mg BID (ASCEMBL trial)
# STAMP mechanism, active against most mutations
#-----------------------------------------------------------------
cat("\n--- Scenario 5: Asciminib 40mg BID (ASCEMBL) ---\n")
scen5 <- mod %>%
param(dose_imt = 0, dose_das = 0, dose_nil = 0, dose_asc = 80) %>%
# 40mg BID = 80mg/day
mrgsim(end = t_end_hr, delta = dt) %>%
as.data.frame() %>%
mutate(scenario = "Asciminib 40mg BID", time_months = to_months(time))
mmr_24mo_asc <- scen5$MMR[scen5$time_months >= 24][1]
cat("MMR @ 24mo (Asciminib):", mmr_24mo_asc, "\n")
#-----------------------------------------------------------------
# SCENARIO 6: Imatinib with T315I resistance emergence
# Initial response, then resistance clone outgrowth
#-----------------------------------------------------------------
cat("\n--- Scenario 6: Imatinib → T315I Resistance Emergence ---\n")
scen6 <- mod %>%
param(dose_imt = 400, dose_das = 0, dose_nil = 0, dose_asc = 0) %>%
# Higher mutation rate for resistance emergence
param(mu = 5e-7) %>%
mrgsim(end = t_end_hr, delta = dt) %>%
as.data.frame() %>%
mutate(scenario = "Imatinib + T315I resistance", time_months = to_months(time))
# When does BCR-ABL IS% start rising (resistance detection)?
first_rise <- which(diff(scen6$BCRPCT) > 0.001 & scen6$BCRPCT[-1] > 1.0)
if (length(first_rise) > 0) {
cat("Resistance detected (IS% rising) at:", round(scen6$time_months[first_rise[1]], 1), "months\n")
}
cat("Resistant fraction @ 5yr:", round(tail(scen6$Resist_frac, 1) * 100, 2), "%\n")
#-----------------------------------------------------------------
# SCENARIO 7: TFR attempt after deep molecular response
# Patient achieves MR4.5 on nilotinib, then stops TKI
#-----------------------------------------------------------------
cat("\n--- Scenario 7: TFR Attempt after Sustained DMR ---\n")
# Phase 1: nilotinib until MR4.5 (simulate 4 years)
t_tfr_start <- 4 * 365.25 * 24 # 4 years in hours
t_tfr_end <- t_end_hr
# First, get state at TFR initiation after 4yr nilotinib
pre_TFR <- mod %>%
param(dose_imt = 0, dose_das = 0, dose_nil = 600, dose_asc = 0) %>%
mrgsim(end = t_tfr_start, delta = dt) %>%
as.data.frame()
cat("BCR-ABL IS% at TFR start (4yr nilotinib):",
round(tail(pre_TFR$BCRPCT, 1), 4), "%\n")
cat("MR4.5 achieved at TFR:", tail(pre_TFR$MR45, 1), "\n")
# Phase 2: TKI stopped — model natural clearance
# Extract end-state for continuation (simplified: reset PK, keep cell states)
end_state <- tail(pre_TFR, 1)
init_tfr <- c(
y0 = end_state$y0,
y1 = end_state$y1,
y2 = end_state$y2,
y3 = end_state$y3,
x0q = end_state$x0q,
x0p = end_state$x0p,
x1 = end_state$x1,
x2 = end_state$x2,
x3 = end_state$x3,
xs = end_state$xs,
xsp = end_state$xsp,
xr = end_state$xr,
xrp = end_state$xrp,
Gut_imt = 0, Cp_imt = 0, Cic_imt = 0,
Cp_das = 0, Cic_das = 0,
Gut_nil = 0, Cp_nil = 0, Cic_nil = 0,
Gut_asc = 0, Cp_asc = 0
)
# Simulate TFR phase (off TKI)
scen7_tfr <- mod %>%
param(dose_imt = 0, dose_das = 0, dose_nil = 0, dose_asc = 0) %>%
init(init_tfr) %>%
mrgsim(end = (6 * 365.25 * 24), delta = dt) %>%
as.data.frame() %>%
mutate(scenario = "TFR (off nilotinib)",
time_months = to_months(time) + 48) # offset by 4yr
# Combine pre- and post-TFR
pre_TFR2 <- pre_TFR %>%
mutate(scenario = "Nilotinib pre-TFR", time_months = to_months(time))
scen7 <- bind_rows(pre_TFR2, scen7_tfr)
# TFR outcome: molecular relapse if IS% > 0.1% within 6 months
tfr_6mo <- scen7_tfr$BCRPCT[scen7_tfr$time_months >= 48 + 6][1]
if (!is.na(tfr_6mo) && tfr_6mo > 0.1) {
cat("TFR outcome: MOLECULAR RELAPSE @ 6mo (IS% =", round(tfr_6mo, 4), "%)\n")
} else {
cat("TFR outcome: IN REMISSION @ 6mo (IS% =", round(tfr_6mo, 4), "%)\n")
}
# ============================================================
# COMBINED PLOT: All TKI scenarios — BCR-ABL IS%
# ============================================================
all_scenarios <- bind_rows(scen1, scen2, scen3, scen4, scen5, scen6) %>%
filter(time_months <= 60) # 5 years
p1 <- ggplot(all_scenarios, aes(x = time_months, y = BCRPCT,
color = scenario, linetype = scenario)) +
geom_line(size = 1.2) +
scale_y_log10(
breaks = c(100, 10, 1, 0.1, 0.01, 0.001),
labels = c("100%", "10%", "1%", "0.1%\n(MMR)", "0.01%\n(MR4)", "0.001%\n(MR4.5)"),
limits = c(0.001, 100)
) +
geom_hline(yintercept = 0.1, linetype = "dashed", color = "blue", alpha = 0.7) +
geom_hline(yintercept = 0.01, linetype = "dashed", color = "purple", alpha = 0.7) +
geom_hline(yintercept = 0.0032, linetype = "dashed", color = "darkgreen", alpha = 0.7) +
geom_vline(xintercept = c(3, 6, 12, 18, 36), linetype = "dotted", alpha = 0.4) +
annotate("text", x = 18, y = 0.12, label = "MMR (0.1%)", color = "blue", size = 3) +
annotate("text", x = 18, y = 0.012, label = "MR4 (0.01%)", color = "purple", size = 3) +
annotate("text", x = 18, y = 0.0035, label = "MR4.5", color = "darkgreen", size = 3) +
scale_color_manual(values = c(
"Untreated" = "#E41A1C",
"Imatinib 400mg/day" = "#377EB8",
"Dasatinib 100mg/day" = "#4DAF4A",
"Nilotinib 300mg BID" = "#984EA3",
"Asciminib 40mg BID" = "#FF7F00",
"Imatinib + T315I resistance" = "#A65628"
)) +
labs(
title = "CML QSP Model — BCR-ABL IS% Over Time by TKI Therapy",
subtitle = "Calibrated to IRIS (imatinib), DASISION (dasatinib), ENESTnd (nilotinib), ASCEMBL (asciminib)",
x = "Time (months)",
y = "BCR-ABL IS% (International Scale)",
color = "Treatment Scenario",
linetype = "Treatment Scenario",
caption = "Dashed lines: ELN 2020 response milestones | Dotted verticals: key assessment timepoints"
) +
theme_bw(base_size = 12) +
theme(legend.position = "right")
# ============================================================
# PLOT 2: Cell population kinetics (Imatinib scenario)
# ============================================================
cells_imt <- scen2 %>%
filter(time_months <= 60) %>%
select(time_months, x0q, x0p, x1, x2, x3, y3, xs, xr) %>%
pivot_longer(cols = -time_months, names_to = "population", values_to = "count") %>%
mutate(population = case_when(
population == "x0q" ~ "LSC Quiescent (CML)",
population == "x0p" ~ "LSC Cycling (CML)",
population == "x1" ~ "LPC1 (CML)",
population == "x2" ~ "LPC2 (CML)",
population == "x3" ~ "Leukemic WBC",
population == "y3" ~ "Normal WBC",
population == "xs" ~ "Sensitive LSC Clone",
population == "xr" ~ "Resistant LSC (T315I)"
))
p2 <- ggplot(cells_imt, aes(x = time_months, y = count + 1,
color = population)) +
geom_line(size = 1.0) +
scale_y_log10() +
labs(
title = "CML Cell Population Kinetics — Imatinib 400mg/day",
x = "Time (months)", y = "Cell Count (log10)",
color = "Cell Population"
) +
theme_bw(base_size = 12) +
theme(legend.position = "right")
# ============================================================
# PLOT 3: PK profiles — Imatinib, Dasatinib, Nilotinib
# ============================================================
pk_data <- bind_rows(
scen2 %>% filter(time_months <= 1) %>%
select(time_months, Cp_imt, Cic_imt) %>%
pivot_longer(-time_months, names_to = "compartment", values_to = "conc") %>%
mutate(drug = "Imatinib"),
scen3 %>% filter(time_months <= 1) %>%
select(time_months, Cp_das, Cic_das) %>%
pivot_longer(-time_months, names_to = "compartment", values_to = "conc") %>%
mutate(drug = "Dasatinib"),
scen4 %>% filter(time_months <= 1) %>%
select(time_months, Cp_nil, Cic_nil) %>%
pivot_longer(-time_months, names_to = "compartment", values_to = "conc") %>%
mutate(drug = "Nilotinib")
)
p3 <- ggplot(pk_data, aes(x = time_months * 30 * 24, y = conc,
color = compartment, linetype = drug)) +
geom_line(size = 1.0) +
scale_y_log10() +
labs(
title = "TKI Pharmacokinetic Profiles (First Month)",
x = "Time (hours)", y = "Concentration (µM, log scale)",
color = "Compartment", linetype = "Drug"
) +
theme_bw(base_size = 12)
# Print summary
cat("\n=== SIMULATION COMPLETE ===\n")
cat("Scenarios simulated:\n")
cat(" 1. Untreated CML (progression model)\n")
cat(" 2. Imatinib 400mg/day (IRIS calibration)\n")
cat(" 3. Dasatinib 100mg/day (DASISION calibration)\n")
cat(" 4. Nilotinib 300mg BID (ENESTnd calibration)\n")
cat(" 5. Asciminib 40mg BID (ASCEMBL calibration)\n")
cat(" 6. Imatinib + T315I resistance emergence\n")
cat(" 7. TFR attempt after sustained DMR on nilotinib\n")
cat("\nPlots created: p1 (BCR-ABL IS%), p2 (cell kinetics), p3 (PK)\n")
cat("Use print(p1), print(p2), print(p3) to display\n")
# Return models and results invisibly
invisible(list(
model = mod,
scenarios = list(
untreated = scen1,
imatinib = scen2,
dasatinib = scen3,
nilotinib = scen4,
asciminib = scen5,
resistance = scen6,
tfr = scen7
),
plots = list(p1 = p1, p2 = p2, p3 = p3)
))