Repository navigation
Expand file tree
/
Copy pathmogad_mrgsolve_model.R
More file actions
743 lines (695 loc) · 39.9 KB
/
Copy pathmogad_mrgsolve_model.R
File metadata and controls
743 lines (695 loc) · 39.9 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
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
## MOG Antibody-Associated Disease (MOGAD) — QSP / PK-PD Model
## ============================================================================
## Disease : MOGAD — a CNS-restricted, antibody-mediated demyelinating disease
## in which class-switched IgG1 against myelin oligodendrocyte
## glycoprotein (MOG) reaches the CNS through a permeabilised
## blood-brain barrier, opsonises the OUTERMOST myelin lamella, and
## destroys myelin by complement-dependent lysis and
## antibody-dependent cellular phagocytosis — while sparing the
## oligodendrocyte cell body and the astrocyte/AQP4 compartment.
## Attacks are optic neuritis, myelitis, ADEM and cortical
## encephalitis (FLAMES); disability accrues step-wise from relapses
## with no secondary progressive phase.
## Authors : Generated by Claude Code Routine (QSP Disease Model Library)
## Engine : mrgsolve (>= 1.5) — mod <- mread_cache("mogad_mrgsolve_model.R")
##
## ---------------------------------------------------------------------------
## THREE STRUCTURAL COMMITMENTS
## ---------------------------------------------------------------------------
## (1) SOURCE. The pathogenic IgG1 is made by a SHORT-LIVED, CD20-NEGATIVE
## plasmablast pool (t1/2 = 1/KOUT_PB = 7.1 d) whose generation is only
## partly fed by CD20+ memory B cells. FRAC_ESC (= 0.60) is the share of
## plasmablast generation that rituximab cannot reach. Everything
## downstream of the plasmablast inherits that short time constant, which
## is why titres fall within days of steroids/PLEX and rebound within
## weeks of stopping them.
##
## (2) TARGET. MYEL_* (myelin) and AXON_* (axon) are separate states, and OL
## (surviving oligodendrocytes) multiplies the remyelination rate. In
## MOGAD, KOLD is small — myelin is destroyed, the oligodendrocyte lives,
## so remyelination proceeds and the recovery ceiling is high. Setting
## AQP4_MODE = 1 multiplies oligodendrocyte death by (1 + AQP4_OLD) and
## astrocytic GFAP release by (1 + AQP4_GFAP), converting the SAME model
## into an AQP4-IgG NMOSD comparator with poor recovery and high sGFAP.
## This one switch is the model's built-in control experiment.
##
## (3) NODE OF ACTION. Relapse hazard is a Hill function (HAZ_H = 3) of the
## FUNCTIONALLY ACTIVE antibody AB_EFF, not of the measured titre. Drugs
## act either at the antibody SOURCE (steroids, MMF/AZA, rituximab, IL-6R
## blockade) or at the CLEARANCE/EFFECTOR node (PLEX, IVIG, FcRn
## inhibition, C5 inhibition). Because baseline titre sits on the Hill
## plateau, partial source suppression buys little, and the simulated
## ordering of therapies does NOT follow the simulated ordering of titre
## reduction — see the calibration notes at the bottom.
##
## Scope : 34 ODE compartments —
## * drug PK/exposure (11): methylprednisolone, oral prednisone
## (gut+central) and a GC effect compartment, IVIG, rituximab,
## IL-6R blocker, antimetabolite, FcRn inhibitor, C5 inhibitor,
## a plasma-exchange driver, and an attack-trigger driver
## * immune source (7): CD20+ memory B cells, plasmablasts, serum
## MOG-IgG, total serum IgG, IL-6, activated T cells, BBB
## permeability
## * CNS effector arm (3): parenchymal MOG-IgG, complement
## activation, activated microglia/macrophages
## * tissue (6): optic-nerve and cord myelin and axons, surviving
## oligodendrocytes, lesion oedema
## * biomarkers (3): serum NfL, serum GFAP, CSF cell count
## * cumulative (4): relapse hazard, prednisone-equivalent dose,
## bone mineral density
## Time : days. Default horizon 365 d (index attack at t = 0, IV
## methylprednisolone from day 4, oral prednisone from day 9,
## maintenance therapy from day 28).
## Units : antibody titre in arbitrary units mapped to a live cell-based
## assay dilution via TITER_REF (AB = 3 -> 1:300); myelin/axon/
## oligodendrocyte states as surviving fractions (0-1); IgG in g/L;
## NfL and GFAP in pg/mL; CSF cells in cells/uL; prednisone in mg
## and cumulative dose in g; hazard in 1/day (ARR = 365 * hazard).
## Scenarios: 17 prebuilt — see MOGAD_simulate_scenarios() at the bottom.
## Verify : `mogad_reference_impl.py` in this directory re-implements these
## same 34 ODEs with the same parameter values in dependency-free
## Python (RK4) and prints every number quoted in README.md and in
## the calibration notes below. Run `python3 mogad_reference_impl.py`.
## Disclaimer: For research / education / hypothesis generation only.
## Not a substitute for clinical judgement or regulatory filings.
## ============================================================================
[PROB]
QSP model of MOG antibody-associated disease (MOGAD) linking a short-lived
CD20-negative plasmablast source of MOG-IgG1, FcRn-dependent IgG turnover,
IL-6-driven blood-brain-barrier opening, complement- and
phagocytosis-mediated destruction of the outer myelin lamella with
oligodendrocyte survival, remyelination-dependent recovery, serum NfL/GFAP
biomarkers, and a titre-threshold relapse hazard — with PK/PD for IV
methylprednisolone, oral prednisone tapers, plasma exchange, acute and
maintenance IVIG, rituximab, IL-6 receptor blockade, antimetabolites, FcRn
inhibition and C5 inhibition, plus an AQP4-IgG NMOSD comparator switch.
[PLUGIN] Rcpp
[PARAM] @annotated
//--- Patient / phenotype ----------------------------------------------------
SEV : 1.00 : Attack drive scaling (0.5 mild .. 1.5 severe)
SITE_ON : 1.00 : Topography weight, optic nerve
SITE_SC : 0.45 : Topography weight, spinal cord
AQP4_MODE : 0 : 1 = AQP4-IgG NMOSD comparator arm (astrocytopathy ON)
AQP4_OLD : 25.0 : Extra oligodendrocyte-death coefficient when AQP4_MODE=1
AQP4_GFAP : 30.0 : Extra astrocytic GFAP release when AQP4_MODE=1
BYST_GFAP : 0.60 : Bystander GFAP release in MOGAD (small by construction)
//--- Glucocorticoids: IV methylprednisolone + oral prednisone taper ---------
KEL_MP : 6.60 : 1/day, methylprednisolone elimination (t1/2 ~2.5 h)
MPPOT : 1.20 : MP -> prednisone-equivalent potency factor
KA_PRED : 8.00 : 1/day, oral prednisone absorption
F_PRED : 0.85 : Oral prednisone bioavailability
KEL_PRED : 5.50 : 1/day, prednisolone elimination (t1/2 ~3 h)
PRED0 : 0 : mg/day, starting oral prednisone dose (0 = none)
TAPER_DAYS : 1.00 : days, length of the linear oral prednisone wean
PRED_START : 0 : day oral prednisone begins
KIN_GC : 0.70 : 1/day, GC effect-compartment input
KOUT_GC : 0.70 : 1/day, GC effect-compartment turnover
GC50 : 4.70 : GC effect units giving half-maximal receptor occupancy
IMAX_GC_PB : 0.85 : Max (reversible) plasmablast suppression by GC
IMAX_GC_IL6: 0.60 : Max IL-6 transcription suppression by GC
IMAX_GC_T : 0.55 : Max T-cell activation suppression by GC
IMAX_GC_BBB: 0.60 : Max fractional increase in BBB resealing rate by GC
IMAX_GC_M : 0.35 : Max myeloid activation suppression by GC
//--- IVIG (acute 2 g/kg; maintenance 1 g/kg q4wk) ---------------------------
KEL_IVIG : 0.040 : 1/day, exogenous IgG elimination (t1/2 ~17 d)
EC50_IVIG : 4.00 : g/L above baseline for half-maximal IVIG effect
FRN_IVIG : 0.45 : Max fractional acceleration of endogenous IgG catabolism
KNEUT : 1.35 : IVIG effector-level neutralisation potency (FcgR blockade, anti-idiotype, complement scavenging)
//--- Rituximab (anti-CD20; SOURCE node) ------------------------------------
KEL_RTX : 0.075 : 1/day, rituximab exposure decay
EC50_RTX : 0.25 : au, half-maximal CD20 occupancy
KDEP_RTX : 0.60 : 1/day, max CD20+ memory B-cell depletion rate
//--- IL-6 receptor blocker (tocilizumab / satralizumab) --------------------
KEL_TCZ : 0.035 : 1/day, IL-6R blocker decay (t1/2 ~20 d)
EC50_TCZ : 0.15 : au, half-maximal IL-6R occupancy
IMAX_TCZ : 0.85 : Max fractional IL-6 receptor blockade
//--- Antimetabolite (mycophenolate / azathioprine) -------------------------
KEL_MMF : 0.35 : 1/day, antimetabolite exposure decay
EC50_MMF : 0.60 : au, half-maximal antiproliferative effect
IMAX_MMF : 0.35 : Max antiproliferative effect
MMF_RATE : 0 : au/day, continuous antimetabolite input (0 = off)
//--- FcRn inhibitor (rozanolixizumab class; CLEARANCE node) ---------------
KEL_FCRN : 0.30 : 1/day, FcRn inhibitor decay
EC50_FCRN : 0.40 : au, half-maximal FcRn occupancy
FRN_INH : 2.30 : Max fold-increase in IgG catabolism at full FcRn blockade
//--- C5 inhibitor (eculizumab class; EFFECTOR node) -----------------------
KEL_C5I : 0.060 : 1/day, C5 inhibitor decay
EC50_C5I : 0.35 : au, half-maximal C5 blockade
IMAX_C5I : 0.95 : Max terminal-pathway (MAC) blockade
IMAX_C5H : 0.35 : Max hazard reduction attributable to MAC blockade
//--- Plasma exchange -------------------------------------------------------
KOUT_PLEX : 4.00 : 1/day, decay of the per-exchange removal driver
KPLEX : 0.90 : 1/day per unit driver; 5 exchanges remove ~2/3 of IgG
//--- Attack trigger (intercurrent infection / index attack) ---------------
KOUT_TRIG : 0.30 : 1/day, trigger decay (t1/2 ~2.3 d)
ATRIG : 1.50 : Trigger -> memory B-cell expansion gain
BTRIG : 6.00 : Trigger -> plasmablast burst gain
CTRIG : 25.0 : Trigger -> IL-6 gain
DTRIG : 5.00 : Trigger -> T-cell activation gain
ETRIG : 70.0 : Trigger -> BBB opening gain
//--- B lineage / antibody (commitment 1) ----------------------------------
KIN_MBC : 0.008 : 1/day, memory B-cell replenishment
KOUT_MBC : 0.008 : 1/day, memory B-cell loss (repopulation t1/2 ~87 d)
KGEN_PB : 0.030435 : 1/day, plasmablast generation rate constant
KOUT_PB : 0.140 : 1/day, plasmablast loss (t1/2 ~7 d — SHORT-LIVED)
FRAC_ESC : 0.60 : CD20-negative share of plasmablast generation (the share rituximab structurally cannot reach)
KSYN_AB : 0.099 : au/day per plasmablast unit, MOG-IgG1 secretion
KEL_AB : 0.033 : 1/day, IgG elimination (t1/2 ~21 d, FcRn-set)
KSYN_IgG : 0.330 : g/L/day, total IgG synthesis
FLOOR_IgG : 0.40 : Fraction of total IgG synthesis independent of CD20+ cells
EMAX_IL6 : 6.00 : Max fold-gain of plasmablast generation by IL-6 signalling
EC50_IL6 : 2.00 : pg/mL IL-6 signal for half-maximal plasmablast support
//--- Cytokines / T cells / barrier ---------------------------------------
KIN_IL6 : 8.00 : pg/mL/day, IL-6 production
KOUT_IL6 : 4.00 : 1/day, IL-6 turnover
KT_IL6 : 0.50 : T-cell drive on IL-6 production
KIN_T : 0.080 : 1/day, T-cell activation input
KOUT_T : 0.080 : 1/day, activated T-cell loss (t1/2 ~8.7 d)
KIN_BBB : 0.00918 : 1/day, BBB opening rate constant
KOUT_BBB : 0.50 : 1/day, BBB resealing
KIL6_BBB : 2.00 : IL-6 drive on BBB opening
EC50_IL6B : 6.00 : pg/mL IL-6 for half-maximal BBB effect
KT_BBB : 1.20 : T-cell drive on BBB opening
//--- CNS transit & effectors ---------------------------------------------
KIN_CNS : 0.050 : 1/day, BBB-gated IgG influx into CNS parenchyma
KOUT_CNS : 0.250 : 1/day, CNS IgG clearance (CSF bulk flow)
KIN_COMP : 1.00 : 1/day, classical-pathway activation per unit opsonisation
KOUT_COMP : 1.00 : 1/day, complement decay
KIN_MPHI : 0.30 : 1/day, myeloid activation input
KOUT_MPHI : 0.30 : 1/day, myeloid deactivation
KC5A : 2.00 : C5a chemotactic amplification of myeloid recruitment
KCDC : 1.50 : Weight of complement-dependent myelin lysis
KADCP : 8.00 : Weight of antibody-dependent cellular phagocytosis
KDIR : 0.50 : Weight of direct (antibody-only) myelin destabilisation
INJMAX : 0.90 : 1/day, maximal demyelination rate
OPS50 : 1.20 : Effector signal giving half-maximal demyelination
HI : 2.50 : Hill coefficient of the demyelination threshold
//--- Tissue (commitment 2) ----------------------------------------------
KREMY : 0.035 : 1/day, remyelination rate at full oligodendrocyte survival
KOLD : 0.080 : Oligodendrocyte death per unit injury (SMALL in MOGAD)
KOLREP : 0.010 : 1/day, OPC-driven oligodendrocyte replacement
KAX : 0.220 : Axonal loss rate above the injury threshold
AXTHR : 0.060 : 1/day injury threshold below which axons are spared
KAX2 : 0.010 : Chronic degeneration rate of persistently denuded axons
DEMYTHR : 0.350 : Demyelination fraction above which chronic loss begins
KIN_ED : 3.00 : Oedema formation per unit injury
KOUT_ED : 0.150 : 1/day, oedema resolution
//--- Biomarkers ----------------------------------------------------------
KIN_NFL : 900.0 : pg/mL per unit axonal-loss rate
KBASE_NFL : 0.640 : pg/mL/day, constitutive NfL release (baseline 8 pg/mL)
KOUT_NFL : 0.080 : 1/day, serum NfL clearance
KIN_GFAP : 110.0 : pg/mL per unit astrocytic injury
KBASE_GFAP : 6.40 : pg/mL/day, constitutive GFAP release (baseline 80 pg/mL)
KOUT_GFAP : 0.080 : 1/day, serum GFAP clearance
KIN_CSF : 25.0 : cells/uL/day per unit barrier+myeloid signal
KOUT_CSF : 0.120 : 1/day, CSF cell clearance
KM_CSF : 1.00 : Myeloid contribution to CSF pleocytosis
//--- Relapse hazard (commitment 3) --------------------------------------
LAM0 : 0.00213 : 1/day, maximal relapse hazard (ARR_max = 0.78/y)
AB50 : 1.80 : au, functionally active antibody at half-maximal hazard
HAZ_H : 3.00 : Hill coefficient — the titre THRESHOLD
KAMP : 6.00 : IL-6 gain on attack amplification competence
EC50_AMP : 1.20 : pg/mL IL-6 signal for half-maximal amplification
POW_BBB : 0.70 : Exponent of the BBB term in the hazard
IL6_BASE : 3.00 : pg/mL, baseline IL-6 (hazard normalisation reference)
BBB_BASE : 0.050 : Baseline BBB permeability index (normalisation reference)
//--- Glucocorticoid toxicity -------------------------------------------
KBMD : 2.70e-4 : 1/day BMD loss at full GC receptor occupancy
KBMDREC : 1.00e-4 : 1/day BMD recovery toward baseline off steroid
//--- Endpoint mapping --------------------------------------------------
WM : 2.00 : logMAR per unit optic-nerve demyelination (reversible)
WA : 3.50 : logMAR per unit axonal loss beyond functional reserve
AXRES : 0.15 : Axonal functional reserve (fraction lost before VA falls)
RNFL0 : 98.0 : um, baseline peripapillary RNFL thickness
TITER_REF : 100.0 : Live cell-based assay dilution per unit AB
[CMT] @annotated
MP_C : IV methylprednisolone, prednisone-equivalent exposure (mg)
PRED_G : Oral prednisone gut depot (mg)
PRED_C : Prednisolone central exposure (mg-equivalent)
GC_EFF : Glucocorticoid effect compartment (au)
IVIG_C : Exogenous IVIG-derived IgG above baseline (g/L)
RTX_C : Rituximab exposure (au)
TCZ_C : IL-6 receptor blocker exposure (au)
MMF_C : Antimetabolite exposure (au)
FCRN_C : FcRn inhibitor exposure (au)
C5I_C : C5 inhibitor exposure (au)
PLEX_D : Plasma-exchange removal driver (au)
TRIG : Attack-trigger driver, e.g. intercurrent infection (au)
MBC : CD20+ MOG-specific memory B cells (relative to baseline)
PB : Short-lived CD20-negative plasmablasts (relative to baseline)
IL6 : Interleukin-6 (pg/mL)
TCELL : Activated MOG-specific T cells (relative to baseline)
AB : Serum MOG-IgG1 titre (au; x TITER_REF = CBA dilution)
TOTIgG : Total endogenous serum IgG (g/L)
BBB : Blood-brain-barrier permeability index (0-1)
AB_CNS : CNS parenchymal MOG-IgG (au)
COMP : Local complement activation / MAC deposition (au)
MPHI : Activated microglia / macrophages (au)
MYEL_ON : Optic-nerve myelin integrity (surviving fraction)
AXON_ON : Optic-nerve axonal integrity (surviving fraction)
MYEL_SC : Spinal-cord myelin integrity (surviving fraction)
AXON_SC : Spinal-cord axonal integrity (surviving fraction)
OL : Surviving oligodendrocyte fraction (sets remyelination ceiling)
EDEMA : Lesion oedema / inflammatory volume (au)
NFL : Serum neurofilament light (pg/mL)
GFAP : Serum glial fibrillary acidic protein (pg/mL)
CSFC : CSF white cell count (cells/uL)
CUMHAZ : Cumulative relapse hazard (integral of hazard; P = 1 - exp(-CUMHAZ))
CUMSTER : Cumulative prednisone-equivalent glucocorticoid dose (g)
BMD : Bone mineral density (fraction of baseline)
[MAIN]
MBC_0 = 1.0;
PB_0 = 1.0;
IL6_0 = IL6_BASE;
TCELL_0 = 1.0;
AB_0 = 3.0;
TOTIgG_0 = 10.0;
BBB_0 = BBB_BASE;
AB_CNS_0 = 0.030;
COMP_0 = 0.030;
MPHI_0 = 0.090;
MYEL_ON_0 = 0.97;
AXON_ON_0 = 1.0;
MYEL_SC_0 = 0.97;
AXON_SC_0 = 1.0;
OL_0 = 1.0;
NFL_0 = 8.0;
GFAP_0 = 80.0;
CSFC_0 = 3.0;
BMD_0 = 1.0;
[ODE]
//== Glucocorticoid PK and receptor occupancy ==============================
// The oral prednisone wean is a zero-order input with a LINEAR taper, so
// TAPER_DAYS is a continuously tunable lever and total dose = PRED0*TAPER/2.
double trel = SOLVERTIME - PRED_START;
double PRED_IN = 0.0;
if (PRED0 > 0.0 && trel >= 0.0 && trel < TAPER_DAYS)
PRED_IN = PRED0 * (1.0 - trel / TAPER_DAYS);
dxdt_MP_C = -KEL_MP * MP_C;
dxdt_PRED_G = PRED_IN - KA_PRED * PRED_G;
dxdt_PRED_C = KA_PRED * PRED_G * F_PRED - KEL_PRED * PRED_C;
double PEQ = PRED_C + MPPOT * MP_C;
dxdt_GC_EFF = KIN_GC * PEQ - KOUT_GC * GC_EFF;
double GC_occ = GC_EFF / (GC50 + GC_EFF);
//== Other drug exposures =================================================
dxdt_IVIG_C = -KEL_IVIG * IVIG_C;
dxdt_RTX_C = -KEL_RTX * RTX_C;
dxdt_TCZ_C = -KEL_TCZ * TCZ_C;
dxdt_MMF_C = MMF_RATE - KEL_MMF * MMF_C;
dxdt_FCRN_C = -KEL_FCRN * FCRN_C;
dxdt_C5I_C = -KEL_C5I * C5I_C;
dxdt_PLEX_D = -KOUT_PLEX * PLEX_D;
dxdt_TRIG = -KOUT_TRIG * TRIG;
double RTX_occ = RTX_C / (EC50_RTX + RTX_C);
double IL6R_block = IMAX_TCZ * TCZ_C / (EC50_TCZ + TCZ_C);
double MMF_eff = IMAX_MMF * MMF_C / (EC50_MMF + MMF_C);
double FCRN_occ = FCRN_C / (EC50_FCRN + FCRN_C);
double C5I_block = IMAX_C5I * C5I_C / (EC50_C5I + C5I_C);
double IVIG_sat = IVIG_C / (EC50_IVIG + IVIG_C);
//== COMMITMENT 1: short-lived CD20-negative plasmablast source ============
double IL6_sig = IL6 * (1.0 - IL6R_block);
double FIL6 = 1.0 + EMAX_IL6 * IL6_sig / (EC50_IL6 + IL6_sig);
double GCsup = IMAX_GC_PB * GC_occ;
dxdt_MBC = KIN_MBC * (1.0 - MMF_eff) * (1.0 + ATRIG * TRIG)
- KOUT_MBC * MBC - KDEP_RTX * RTX_occ * MBC;
// FRAC_ESC of plasmablast generation is CD20-independent: rituximab drives
// MBC to ~0 but can never take the bracket below FRAC_ESC.
double gen = (1.0 - FRAC_ESC) * MBC + FRAC_ESC;
dxdt_PB = KGEN_PB * gen * FIL6 * (1.0 - GCsup) * (1.0 - 0.5 * MMF_eff)
* (1.0 + BTRIG * TRIG) - KOUT_PB * PB;
// FcRn-dependent catabolism: IVIG saturates FcRn, FcRn inhibitors block it.
double FcRn_acc = FRN_IVIG * IVIG_sat + FRN_INH * FCRN_occ;
dxdt_AB = KSYN_AB * PB - KEL_AB * AB * (1.0 + FcRn_acc) - KPLEX * PLEX_D * AB;
dxdt_TOTIgG = KSYN_IgG * (FLOOR_IgG + (1.0 - FLOOR_IgG) * MBC)
- KEL_AB * TOTIgG * (1.0 + FcRn_acc) - KPLEX * PLEX_D * TOTIgG;
//== Cytokines, T cells, barrier ==========================================
dxdt_IL6 = KIN_IL6 * (1.0 + CTRIG * TRIG + KT_IL6 * TCELL)
* (1.0 - IMAX_GC_IL6 * GC_occ) - KOUT_IL6 * IL6;
dxdt_TCELL = KIN_T * (1.0 + DTRIG * TRIG) * (1.0 - IMAX_GC_T * GC_occ)
* (1.0 - MMF_eff) - KOUT_T * TCELL;
double bbb_drive = 1.0 + ETRIG * TRIG
+ KIL6_BBB * IL6_sig / (EC50_IL6B + IL6_sig)
+ KT_BBB * TCELL;
dxdt_BBB = KIN_BBB * bbb_drive * (1.0 - BBB)
- KOUT_BBB * BBB * (1.0 + IMAX_GC_BBB * GC_occ);
//== CNS transit and effector arm =========================================
// AB_EFF is the FUNCTIONALLY ACTIVE antibody. IVIG lowers it far more than it
// lowers the measured titre (FcgR blockade, anti-idiotype, C scavenging) —
// this is why titre is a poor surrogate for protection (commitment 3).
double AB_eff = AB / (1.0 + KNEUT * IVIG_sat);
dxdt_AB_CNS = KIN_CNS * BBB * AB_eff - KOUT_CNS * AB_CNS;
double OPS = AB_CNS;
dxdt_COMP = KIN_COMP * OPS * (1.0 - C5I_block) - KOUT_COMP * COMP;
dxdt_MPHI = KIN_MPHI * (OPS + KC5A * COMP) * (1.0 - IMAX_GC_M * GC_occ)
- KOUT_MPHI * MPHI;
double OPS_eff = KCDC * COMP + KADCP * MPHI * OPS + KDIR * OPS;
double ops_h = pow(fmax(OPS_eff, 0.0), HI);
double INJ = SEV * INJMAX * ops_h / (pow(OPS50, HI) + ops_h);
double INJ_ON = SITE_ON * INJ;
double INJ_SC = SITE_SC * INJ;
//== COMMITMENT 2: myelin is lost, the oligodendrocyte survives ============
// Remyelination is proportional to OL, so OL sets the recovery CEILING.
dxdt_MYEL_ON = -INJ_ON * MYEL_ON + KREMY * OL * (1.0 - MYEL_ON);
dxdt_MYEL_SC = -INJ_SC * MYEL_SC + KREMY * OL * (1.0 - MYEL_SC);
double ax_on = KAX * fmax(INJ_ON - AXTHR, 0.0) * AXON_ON
+ KAX2 * fmax(1.0 - MYEL_ON - DEMYTHR, 0.0) * AXON_ON;
double ax_sc = KAX * fmax(INJ_SC - AXTHR, 0.0) * AXON_SC
+ KAX2 * fmax(1.0 - MYEL_SC - DEMYTHR, 0.0) * AXON_SC;
dxdt_AXON_ON = -ax_on;
dxdt_AXON_SC = -ax_sc;
dxdt_OL = -KOLD * (1.0 + AQP4_MODE * AQP4_OLD) * INJ * OL + KOLREP * (1.0 - OL);
dxdt_EDEMA = KIN_ED * INJ - KOUT_ED * EDEMA * (1.0 + 2.0 * GC_occ);
//== Biomarkers ===========================================================
dxdt_NFL = KIN_NFL * (ax_on + ax_sc) + KBASE_NFL - KOUT_NFL * NFL;
double ASTRO_INJ = (BYST_GFAP + AQP4_MODE * AQP4_GFAP) * INJ;
dxdt_GFAP = KIN_GFAP * ASTRO_INJ + KBASE_GFAP - KOUT_GFAP * GFAP;
dxdt_CSFC = KIN_CSF * (BBB + KM_CSF * MPHI) - KOUT_CSF * CSFC;
//== COMMITMENT 3: relapse hazard is a THRESHOLD function of AB_EFF ========
double ab_h = pow(fmax(AB_eff, 0.0), HAZ_H);
double Hill_AB = ab_h / (pow(AB50, HAZ_H) + ab_h);
double amp_il6 = (1.0 + KAMP * IL6_sig / (EC50_AMP + IL6_sig))
/ (1.0 + KAMP * IL6_BASE / (EC50_AMP + IL6_BASE));
double amp_bbb = pow(fmax(BBB, 1e-9) / BBB_BASE, POW_BBB);
double AMP = amp_il6 * amp_bbb * (1.0 - IMAX_C5H * C5I_block);
dxdt_CUMHAZ = LAM0 * Hill_AB * AMP;
//== Glucocorticoid toxicity ==============================================
// Elimination flux integrates to total administered dose.
dxdt_CUMSTER = (KEL_PRED * PRED_C + KEL_MP * MP_C) / 1000.0;
dxdt_BMD = -KBMD * GC_occ + KBMDREC * (1.0 - BMD);
[TABLE]
double GC_occ2 = GC_EFF / (GC50 + GC_EFF);
double IL6R_blk2 = IMAX_TCZ * TCZ_C / (EC50_TCZ + TCZ_C);
double IVIG_sat2 = IVIG_C / (EC50_IVIG + IVIG_C);
double C5I_blk2 = IMAX_C5I * C5I_C / (EC50_C5I + C5I_C);
double RTX_occ2 = RTX_C / (EC50_RTX + RTX_C);
double FCRN_occ2 = FCRN_C / (EC50_FCRN + FCRN_C);
double AB_eff2 = AB / (1.0 + KNEUT * IVIG_sat2);
double IL6_sig2 = IL6 * (1.0 - IL6R_blk2);
double ab_h2 = pow(fmax(AB_eff2, 0.0), HAZ_H);
double Hill_AB2 = ab_h2 / (pow(AB50, HAZ_H) + ab_h2);
double amp_il62 = (1.0 + KAMP * IL6_sig2 / (EC50_AMP + IL6_sig2))
/ (1.0 + KAMP * IL6_BASE / (EC50_AMP + IL6_BASE));
double amp_bbb2 = pow(fmax(BBB, 1e-9) / BBB_BASE, POW_BBB);
double HAZ = LAM0 * Hill_AB2 * amp_il62 * amp_bbb2
* (1.0 - IMAX_C5H * C5I_blk2);
// Deficit = REVERSIBLE conduction block from demyelination + IRREVERSIBLE
// axonal loss beyond a functional reserve. The reserve term is what allows
// marked RNFL thinning together with good visual recovery — the MOGAD
// structure-function dissociation.
double logMAR = WM * (1.0 - MYEL_ON)
+ WA * fmax((1.0 - AXON_ON) - AXRES, 0.0);
if (logMAR > 3.0) logMAR = 3.0;
double edss = 4.5 * (1.0 - MYEL_SC)
+ 9.0 * fmax((1.0 - AXON_SC) - 0.05, 0.0)
+ 0.48 * fmin(logMAR, 2.0);
if (edss > 9.5) edss = 9.5;
capture MOG_IgG_titer = TITER_REF * AB;
capture MOG_IgG_active = AB_eff2;
capture GC_occupancy = GC_occ2;
capture CD20_occupancy = RTX_occ2;
capture IL6R_blockade = IL6R_blk2;
capture FcRn_occupancy = FCRN_occ2;
capture C5_blockade = C5I_blk2;
capture IL6_measured = IL6 * (1.0 + 3.0 * IL6R_blk2);
capture BBB_permeability = BBB;
capture CNS_MOG_IgG = AB_CNS;
capture Complement_MAC = COMP;
capture Myeloid_activation = MPHI;
capture Myelin_optic_nerve = MYEL_ON;
capture Axon_optic_nerve = AXON_ON;
capture Myelin_cord = MYEL_SC;
capture Axon_cord = AXON_SC;
capture Oligodendrocytes = OL;
capture Lesion_edema = EDEMA;
capture Serum_NfL = NFL;
capture Serum_GFAP = GFAP;
capture CSF_cells = CSFC;
capture RNFL_thickness = RNFL0 * AXON_ON;
capture Total_IgG = TOTIgG + IVIG_C;
capture VA_logMAR = logMAR;
capture VA_snellen_denom = 20.0 * pow(10.0, logMAR);
capture EDSS = edss;
capture Relapse_hazard = HAZ;
capture ARR = HAZ * 365.0;
capture Relapse_prob = 1.0 - exp(-CUMHAZ);
capture Steroid_dose_g = CUMSTER;
capture BMD_fraction = BMD;
capture Plasmablasts = PB;
capture Memory_B_cells = MBC;
## Captured outputs (declared with `capture` in [TABLE] above):
## MOG_IgG_titer MOG_IgG_active GC_occupancy CD20_occupancy IL6R_blockade
## FcRn_occupancy C5_blockade IL6_measured BBB_permeability CNS_MOG_IgG
## Complement_MAC Myeloid_activation Myelin_optic_nerve Axon_optic_nerve
## Myelin_cord Axon_cord Oligodendrocytes Lesion_edema Serum_NfL Serum_GFAP
## CSF_cells RNFL_thickness Total_IgG VA_logMAR VA_snellen_denom EDSS
## Relapse_hazard ARR Relapse_prob Steroid_dose_g BMD_fraction Plasmablasts
## Memory_B_cells
## ============================================================================
## R-side scenario helpers (read by users, not by mrgsolve)
## ============================================================================
##
## library(mrgsolve); library(dplyr); library(tidyr); library(ggplot2)
## mod <- mread_cache("mogad_mrgsolve_model.R")
##
## # ---------- Dosing conventions ----------
## # MP_C : 1 g IV methylprednisolone = amt 1250 (prednisone-equivalent mg)
## # IVIG_C : 1 g/kg IVIG = amt 12 (g/L increment)
## # PLEX_D : one plasma exchange = amt 1
## # RTX_C / TCZ_C / FCRN_C / C5I_C : amt 1 per administration (au)
## # TRIG : one attack trigger = amt 1
## # Oral prednisone is NOT an event: set PRED0 (mg/d), TAPER_DAYS and
## # PRED_START, and the linear wean is generated inside [ODE].
## # Index attack = TRIG amt 1 at t = 0; IVMP from day 4 (treatment delay);
## # oral prednisone from day 9; maintenance therapy from day 28.
##
## ev_attack <- function() ev(amt = 1, cmt = "TRIG", time = 0)
## ev_ivmp <- function(d0 = 4, n = 5)
## ev(amt = 1250, cmt = "MP_C", time = d0, ii = 1, addl = n - 1)
## ev_plex <- function(d0 = 6, n = 5, spacing = 2)
## ev(amt = 1, cmt = "PLEX_D", time = d0, ii = spacing, addl = n - 1)
## ev_ivig_a <- function(d0 = 5)
## ev(amt = 12, cmt = "IVIG_C", time = d0, ii = 1, addl = 1)
## ev_ivig_m <- function(d0 = 28, ii = 28, gkg = 1, end = 365)
## ev(amt = 12 * gkg, cmt = "IVIG_C", time = d0, ii = ii,
## addl = max(floor((end - d0) / ii), 0))
## ev_rtx <- function(d0 = 28, end = 365)
## c(ev(amt = 1, cmt = "RTX_C", time = d0, ii = 14, addl = 1),
## ev(amt = 1, cmt = "RTX_C", time = d0 + 182, ii = 182,
## addl = max(floor((end - d0 - 182) / 182), 0)))
## ev_il6rb <- function(d0 = 28, ii = 28, end = 365)
## ev(amt = 1, cmt = "TCZ_C", time = d0, ii = ii,
## addl = max(floor((end - d0) / ii), 0))
## ev_fcrn <- function(d0 = 28, ii = 7, end = 365)
## ev(amt = 1, cmt = "FCRN_C", time = d0, ii = ii,
## addl = max(floor((end - d0) / ii), 0))
## ev_c5i <- function(d0 = 28, ii = 14, end = 365)
## ev(amt = 1, cmt = "C5I_C", time = d0, ii = ii,
## addl = max(floor((end - d0) / ii), 0))
##
## # ---------- 17 prebuilt scenarios ----------
## MOGAD_simulate_scenarios <- function(end_d = 365) {
##
## # Standard acute treatment + guideline-like 90-day oral wean
## ST <- list(PRED0 = 60, TAPER_DAYS = 90, PRED_START = 9)
## A <- ev_attack(); M <- ev_ivmp(4, 5)
##
## sc <- list(
## "S1 Untreated index attack" =
## list(par = list(), ev = A),
## "S2 IVMP x5 + SHORT taper (60 mg over 28 d; 840 mg total)" =
## list(par = list(PRED0 = 60, TAPER_DAYS = 28, PRED_START = 9), ev = c(A, M)),
## "S3 IVMP x5 + LONG taper (18.7 mg over 90 d; 840 mg total)" =
## list(par = list(PRED0 = 18.67, TAPER_DAYS = 90, PRED_START = 9), ev = c(A, M)),
## "S4 IVMP x5 + guideline taper (60 mg over 90 d; 2700 mg)" =
## list(par = ST, ev = c(A, M)),
## "S5 IVMP + PLEX x5 (severe attack)" =
## list(par = c(ST, list(SEV = 1.35)), ev = c(A, M, ev_plex(6, 5, 2))),
## "S6 IVMP + acute IVIG 2 g/kg" =
## list(par = ST, ev = c(A, M, ev_ivig_a(5))),
## "S7 Maintenance IVIG 1 g/kg q4wk" =
## list(par = ST, ev = c(A, M, ev_ivig_m(28, 28, 1.0, end_d))),
## "S8 Maintenance IVIG 0.4 g/kg q8wk (under-dosed)" =
## list(par = ST, ev = c(A, M, ev_ivig_m(28, 56, 0.4, end_d))),
## "S9 Rituximab (1 g x2, then q6mo)" =
## list(par = ST, ev = c(A, M, ev_rtx(28, end_d))),
## "S10 IL-6R blockade (tocilizumab / satralizumab q4wk)" =
## list(par = ST, ev = c(A, M, ev_il6rb(28, 28, end_d))),
## "S11 Mycophenolate / azathioprine (continuous)" =
## list(par = c(ST, list(MMF_RATE = 0.42)), ev = c(A, M)),
## "S12 FcRn inhibitor (weekly SC)" =
## list(par = ST, ev = c(A, M, ev_fcrn(28, 7, end_d))),
## "S13 C5 inhibitor (q2wk)" =
## list(par = ST, ev = c(A, M, ev_c5i(28, 14, end_d))),
## "S14 Maintenance IVIG + IL-6R blockade (combination)" =
## list(par = ST, ev = c(A, M, ev_ivig_m(28, 28, 1.0, end_d),
## ev_il6rb(28, 28, end_d))),
## "S15 AQP4-IgG NMOSD comparator (same attack, same IVMP)" =
## list(par = c(ST, list(AQP4_MODE = 1)), ev = c(A, M)),
## "S16 Intercurrent infection at day 200 ON maintenance IVIG" =
## list(par = ST, ev = c(A, M, ev_ivig_m(28, 28, 1.0, end_d),
## ev(amt = 1, cmt = "TRIG", time = 200))),
## "S17 Intercurrent infection at day 200, NO maintenance" =
## list(par = ST, ev = c(A, M, ev(amt = 1, cmt = "TRIG", time = 200)))
## )
##
## out <- lapply(names(sc), function(nm) {
## s <- sc[[nm]]
## mod %>%
## param(s$par) %>%
## ev(s$ev) %>%
## mrgsim(end = end_d, delta = 0.5) %>%
## as_tibble() %>%
## mutate(scenario = nm)
## })
## bind_rows(out)
## }
##
## # ---------- Taper-duration sweep at matched cumulative dose ----------
## MOGAD_taper_sweep <- function(total_mg = 840, end_d = 365) {
## tds <- c(14, 28, 45, 60, 90, 120, 180)
## bind_rows(lapply(tds, function(td) {
## mod %>%
## param(PRED0 = 2 * total_mg / td, TAPER_DAYS = td, PRED_START = 9) %>%
## ev(c(ev_attack(), ev_ivmp(4, 5))) %>%
## mrgsim(end = end_d, delta = 1) %>%
## as_tibble() %>%
## mutate(taper_days = td)
## }))
## }
##
## ============================================================================
## CALIBRATION NOTES — what was fitted, to what, and what the model MISSES
## ============================================================================
## Every number below is printed by `python3 mogad_reference_impl.py`.
##
## ---- Anchors used (all from published cohorts) ----
## * Vilaseca 2026 JAMA Neurol (PMID 42440328; Americas MOGAD Treatment Group,
## n = 116 on IL-6R blockade vs n = 59 on IVIG): relapsing MOGAD ARR before
## preventive therapy 0.64 (95% CI 0.58-0.70); ARR on IL-6R blockade 0.09
## (0.06-0.14); ARR on IVIG 0.22 (0.15-0.32); after IPTW, IL-6R blockade was
## NOT significantly different from IVIG >= 1 g/kg q4wk (HR 2.0, 0.8-4.5) but
## was better than lower-dose IVIG (HR 4.5, 2.0-9.8).
## -> LAM0 was set so that the untreated inter-attack ARR is 0.64 exactly.
## * Chen 2022 JAMA Neurol (PMID 35377395): maintenance IVIG; relapse in 5/29
## (17%) at >= 1 g/kg q4wk vs 15/30 (50%) at lower/less frequent dosing
## (HR 3.31). -> reproduced qualitatively: model ARR 0.18 at 1 g/kg q4wk
## vs 0.37 at 0.4 g/kg q8wk (scenarios S7 vs S8).
## * Chen 2020 Neurology (PMID 32554760) and Thakolwiboon 2021 meta-analysis
## (PMID 34634625): on-treatment ARR ~0.59-0.63 (rituximab), ~0.67-0.84
## (mycophenolate), ~0.2-0.29 (azathioprine), ~0.08 (maintenance IVIG);
## relapse still occurs in ~49-74% on the antimetabolites and rituximab.
## * Chen 2018 Am J Ophthalmol (PMID 30055153): MOG-ON — VA at nadir of the
## worst attack averaged count-fingers, average FINAL VA 20/30, only 6%
## ended <= 20/200; disc oedema and pain each 86%; perineural enhancement
## 50%; longitudinally extensive 80%; persistent seropositivity 61/62.
## * Dubey 2019 JAMA Neurol (PMID 30575890): MOG myelitis — ~1/3 wheelchair-
## dependent at nadir yet recovering better than AQP4-IgG myelitis;
## CSF oligoclonal bands in only 1/38.
## * Stiebel-Kalish 2017 PLoS One (PMID 28125740): RNFL better preserved after
## MOG-IgG than AQP4-IgG optic neuritis.
## * Marignier 2025 JNNP (PMID 39939136) / Mariotto 2019 JNNP (PMID 30952681):
## serum NfL rises with attack severity in MOGAD while serum GFAP stays low,
## in contrast to AQP4-IgG NMOSD.
## * Lledo-Garcia 2022 CPT PSP (PMID 34735735): anti-FcRn (rozanolixizumab)
## PK/PD — total IgG reductions of roughly 50-70% at clinical doses; FRN_INH
## was set to give ~53% total-IgG reduction at steady state (10.0 -> 4.7 g/L).
## * Remlinger 2022 NNN (PMID 35027475): anti-FcRn treatment ameliorates
## MOG-IgG-associated EAE — the mechanistic basis for the FcRn arm.
##
## ---- What the model reproduces (365-day horizon, index optic neuritis) ----
## Therapy model ARR observed ARR (source)
## none (steroid wean only) 0.64 0.64 (0.58-0.70) Vilaseca 2026
## maintenance IVIG 1 g/kg 0.18 0.22 (0.15-0.32) Vilaseca 2026
## IL-6R blockade 0.23 0.09 (0.06-0.14) Vilaseca 2026
## rituximab 0.48 0.59-0.63 Chen 2020/meta
## mycophenolate 0.49 0.67-0.84 Chen 2020/meta
## Attack course (S4): VA nadir logMAR 1.60 (~20/800) -> 0.37 at 1 year
## (~20/47); RNFL 98 -> 75 um; serum NfL 8 -> peak 160 pg/mL; serum GFAP
## peak 150 pg/mL; CSF peak 72 cells/uL; surviving oligodendrocytes 0.88.
## AQP4 comparator (S15, AQP4_MODE = 1, identical attack and steroids):
## oligodendrocytes fall to 0.07, VA settles at logMAR 0.79 (~20/122)
## instead of 0.37, RNFL 63 instead of 75 um, and serum GFAP peaks at
## 3661 instead of 150 pg/mL. The MOGAD-vs-NMOSD recovery gap and the
## sGFAP discrimination both fall out of the single KOLD/AQP4_OLD switch.
##
## ---- The result the structure forces (commitment 3) ----
## Titre reduction does NOT predict relapse protection:
## therapy titre vs none ARR vs none
## maintenance IVIG 0.76 0.28
## rituximab 0.70 0.75
## IL-6R blockade 0.59 0.35
## mycophenolate 0.79 0.77
## FcRn inhibitor 0.47 0.39
## C5 inhibitor 1.00 0.75
## Rituximab lowers the measured titre MORE than IVIG does and yet prevents
## far fewer relapses; the C5 inhibitor lowers the titre not at all and still
## helps. Under a threshold hazard (HAZ_H = 3) with the baseline titre on the
## Hill plateau, what matters is whether an intervention takes the
## FUNCTIONALLY ACTIVE antibody past the threshold — so effector- and
## clearance-node drugs outperform source-node drugs at equal titre effect.
## Falsifiable consequence: serial MOG-IgG titres should be a poor surrogate
## for treatment benefit, and IVIG's benefit should survive in patients whose
## titres do not fall. Both are testable in existing cohorts.
##
## ---- A NEGATIVE result worth stating ----
## At matched cumulative dose (840 mg prednisone-equivalent), stretching the
## oral wean from 14 to 180 days moves the 1-year relapse probability only
## from 0.500 to 0.479. The model therefore does NOT support taper duration
## as an important lever in itself: the hazard reduction bought is roughly
## proportional to steroid spent, and a longer wean mostly defers the rebound
## rather than preventing it. Raising the total dose helps a little more
## (2700 mg over 90 d -> 0.46) at a real BMD cost. The mechanistic reading is
## that oral steroids buy time in proportion to exposure and cannot deliver
## remission, which is the argument for starting a steroid-sparing agent
## early rather than for tapering slowly.
##
## ---- Known misses, stated up front ----
## 1. IL-6R BLOCKADE IS UNDER-PREDICTED (0.23 vs observed 0.09). The model
## reaches an IVIG-like effect but not the reported point estimate. Two
## readings, and the model does not distinguish them: (a) the observed
## incidence-rate ratio of 0.08 is a pre-versus-on-treatment comparison in
## a retrospective cohort and is inflated by regression to the mean, and
## the trustworthy contemporaneous comparison (vs IVIG >= 1 g/kg) was NOT
## significant — in which case the model is right; or (b) IL-6R blockade
## has an effect on plasmablast survival or CNS-compartment inflammation
## larger than EMAX_IL6/KAMP allow — in which case the model is wrong.
## Distinguishing them needs the randomised data.
## 2. AZATHIOPRINE vs MYCOPHENOLATE. Both are one antimetabolite arm here,
## so the model cannot reproduce the reported ordering (azathioprine ARR
## ~0.2-0.29 better than mycophenolate ~0.67-0.84). That ordering is
## mechanistically unexplained and may be confounding/channelling in
## retrospective series; the model takes no position.
## 3. FRAC_ESC = 0.60 IS FITTED, NOT MEASURED. It is the single parameter that
## makes rituximab weaker than IVIG, and the whole source-vs-clearance
## conclusion is contingent on it. A competing explanation — that MOGAD
## relapses on rituximab cluster in B-cell repopulation windows rather than
## arising from CD20-negative escape — would be captured by a low FRAC_ESC
## plus under-dosed redosing, and would predict that strict CD19-guided
## redosing rescues rituximab. That is the experiment that separates them.
## 4. ATTACK NADIR IS SOMEWHAT UNDER-PREDICTED. Treated nadir is logMAR 1.60
## (~20/800) against a reported average nadir of count-fingers, and final
## VA is 20/47 against a reported average 20/30. The model is therefore
## pessimistic about recovery and optimistic about nadir; both would
## improve with an explicit conduction-block state separate from myelin
## content, which is not implemented.
## 5. RELAPSES ARE A HAZARD, NOT EVENTS. CUMHAZ gives the probability of
## relapse but no relapse actually fires, so simulated disability never
## accrues step-wise over multiple attacks. Stochastic relapse firing
## (drawing attack times from the hazard and injecting TRIG) is the natural
## extension and is what would be needed for disability trajectories.
## 6. NO PAEDIATRIC ADEM PHYSIOLOGY. Phenotype is set by SITE_ON/SITE_SC
## weights only; the encephalopathy, seizure and cognitive endpoints in the
## mechanistic map have no ODE counterpart, and paediatric ADEM's much
## higher rate of transient seropositivity is not modelled.
## 7. TITRE-TO-ASSAY MAPPING IS LINEAR (TITER_REF). Real live cell-based assay
## dilutions are ordinal and step in twofold increments; comparisons of
## absolute simulated titres to reported dilutions should be read as
## order-of-magnitude only.
##
## # ---------- Quick sanity plots ----------
## # df <- MOGAD_simulate_scenarios(365)
## # ggplot(df, aes(time, VA_logMAR, colour = scenario)) + geom_line()
## # ggplot(df, aes(time, MOG_IgG_titer, colour = scenario)) + geom_line()
## # ggplot(df, aes(time, ARR, colour = scenario)) + geom_line()
## # ggplot(df, aes(time, Oligodendrocytes, colour = scenario)) + geom_line()
## # ggplot(df, aes(time, Serum_GFAP, colour = scenario)) + geom_line()