Repository navigation
Expand file tree
/
Copy pathmsa_mrgsolve_model.R
More file actions
1520 lines (1398 loc) · 73.4 KB
/
Copy pathmsa_mrgsolve_model.R
File metadata and controls
1520 lines (1398 loc) · 73.4 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
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
## =====================================================================
## Multiple System Atrophy (MSA) — QSP model for mrgsolve
## ---------------------------------------------------------------------
## ORGANISING IDEA
## ---------------
## Upright blood pressure is written as a PRODUCT of independently
## breakable gains, not as a rate:
##
## dMAP(upright) = [baroreceptor afferent] (intact in MSA)
## x [central integration NTS/RVLM] (G_CENT: dying)
## x [preganglionic IML output] (G_CENT: dying)
## x [postganglionic releasable NE] (POSTG: SPARED)
## x [vascular alpha1 density] (A1R: UP-regulated)
## x [effective circulating volume] (ECF: eroded at night)
## x [venous capacitance] (CAPREL)
##
## MSA breaks the CENTRAL terms and spares the POSTGANGLIONIC one.
## Pure autonomic failure (PAF) and Parkinson disease (PD) do the reverse.
## Every phenotype in this file is produced by changing only the VULN_*
## vulnerability parameters — never by changing an equation. Four clinical
## facts then EMERGE rather than being coded:
##
## (1) DRUG SELECTIVITY IS A MAP OF *WHERE* THE PRODUCT IS BROKEN.
## Atomoxetine enters only as a multiplier on NET clearance of NEs, so
## its pressor effect is proportional to NEREL = SNA x NEVES x POSTG.
## In MSA (POSTG high, SNA low-but-present) it works; in the PAF
## phenotype (POSTG -> 0) the same dose does essentially nothing.
## Midodrine's active metabolite is added to the agonist term DIRECTLY,
## downstream of every lesion, so it works in both. Nothing in the
## code says "atomoxetine is MSA-selective".
##
## (2) CARBIDOPA ANTAGONISES DROXIDOPA. Droxidopa becomes NE only through
## DROXNE = KAADC*CDRX*POSTG*(1-CARBI). Carbidopa, given for the
## parkinsonism, occupies the SAME peripheral AADC term. One shared
## enzyme, written once, produces the interaction.
##
## (3) THE PRESSOR CEILING IS SET AT NIGHT, NOT BY DAYTIME POTENCY.
## Supine hypertension -> pressure natriuresis (UNAV, UV both rise
## exponentially in MAP) -> overnight loss of NAB/ECF -> smaller Vc on
## standing next morning -> WORSE orthostatic hypotension. Because the
## loop is closed, bedtime desmopressin and head-up-tilt sleeping raise
## DAYTIME standing BP with no daytime pressor at all, and a long-acting
## pressor at bedtime makes the next morning worse.
##
## (4) LEVODOPA FAILURE IN MSA-P IS POSTSYNAPTIC. Striatal output is
## STRIAT = G_POST * DA/(EC50DA+DA) with G_POST = NMSN^GEXP_MSN.
## Levodopa can only raise DA. In the PD phenotype NMSN stays at 1 and
## levodopa nearly normalises STRIAT; in MSA-P NMSN falls and the SAME
## brain exposure buys a shrinking benefit that wanes over 1-2 years.
##
## TWO TIMESCALES, ONE MODEL
## -------------------------
## Time unit is HOURS. Neurodegeneration has time constants of years;
## the baroreflex has time constants of seconds (TAU_MAP = 0.002 h ~ 7 s).
## Rather than fight the stiffness in one long run, the driver below uses a
## two-stage architecture:
## msa_history() — POSMODE=2 (circadian), 12 years, 12-h output
## msa_state_at() — freeze the state vector at a chosen disease year
## msa_tilt() — restart from that state, POSMODE=3, 40-min tilt at
## 30-s resolution -> DeltaSBP(3 min), DeltaHR/DeltaSBP
## msa_day() — restart from that state, 48-h circadian run with
## real dosing -> supine HTN, nocturnal natriuresis
## The tilt test therefore becomes a READOUT OF DISEASE STAGE computed from
## the same 61 equations that generate the 12-year trajectory.
##
## UNITS
## -----
## time h · volumes L · MAP/SBP mmHg · HR bpm · sodium mmol · plasma drug
## concentrations ng/mL (levodopa/brain levodopa ug/mL) · antibody mg/L ·
## NES/NEVES/SNA/TPR/A1R dimensionless (1 = healthy supine reference) ·
## neuronal states are SURVIVING FRACTIONS (1 = intact).
##
## CALIBRATION TARGETS (see msa_references.md for the source of each)
## ------------------------------------------------------------------
## * nOH definition: standing SBP fall >= 30 mmHg (or DBP >= 15) at 3 min
## * neurogenic signature: DeltaHR/DeltaSBP < 0.5 bpm/mmHg
## * supine plasma NE normal (~250 pg/mL) with a blunted upright increment
## (< ~60%, vs > 100% in health) — the MSA "central lesion" fingerprint
## * UMSARS-II progresses ~5-8 points/year (natural-history cohorts)
## * median survival 6-10 years from motor onset
## * midodrine 10 mg: +15 to +25 mmHg standing SBP, peak ~1 h, ~4 h duration
## * droxidopa 300 mg t.i.d.: ~ +8 to +12 mmHg standing SBP, OHSA item-1
## improvement of ~0.9-1.0 units
## * atomoxetine 18 mg: +25 to +50 mmHg seated SBP in MSA, ~nil in PAF
## * fludrocortisone 0.1 mg: volume-mediated, days to plateau
## * supine hypertension in roughly half of treated patients
## * ~30% of MSA-P show a partial, waning levodopa response; no dyskinesia
## * verdiperstat phase 3 (M-STAR) NEGATIVE — the model must reproduce a
## null, not a benefit, when the intervention starts at diagnosis
##
## RUN
## ---
## Rscript msa_mrgsolve_model.R # builds, self-tests, 19 scenarios
##
## Educational / research model. NOT for clinical use.
## =====================================================================
suppressPackageStartupMessages({
library(mrgsolve)
})
msa_code <- '
$PROB
# Multiple System Atrophy — 61-ODE QSP model (time unit: hours)
$GLOBAL
#define SIGM(x, k) (1.0/(1.0 + exp(-(k)*(x))))
// smooth Gaussian meal bumps at 08:00, 13:00, 19:00
#define MEALB(t) (exp(-pow(((t)-8.0)/1.1,2)) + exp(-pow(((t)-13.0)/1.1,2)) + exp(-pow(((t)-19.0)/1.1,2)))
$PARAM @annotated
// ---------- protocol / posture ----------
POSMODE : 2 : 0 supine, 1 upright, 2 circadian, 3 tilt protocol
UPFRAC : 0.55 : fraction of waking hours spent upright (circadian mode)
HUTF : 0.0 : night-time orthostatic load from head-up tilt sleeping (0-0.3)
WAKE_H : 7 : wake time (h of day)
SLEEP_H : 22 : bedtime (h of day)
TILT_T : 0.25 : time of head-up tilt in tilt mode (h)
MEALAMP : 1.0 : post-prandial splanchnic shunt amplitude (0 = fasting)
BINDER : 0 : abdominal binder / compression garment (0-1)
COUNTERM: 0 : physical counter-manoeuvres + slow rising (0-1)
// ---------- phenotype: regional vulnerability (this IS the phenotype) ----
KND : 0.13 : neurodegeneration drive scale (per year at unit pathology)
VULN_IML : 1.00 : IML preganglionic sympathetic vulnerability
VULN_MED : 0.85 : NTS/CVLM/RVLM medullary cardiovascular vulnerability
VULN_VAG : 0.45 : cardiovagal (nucleus ambiguus) vulnerability
VULN_SN : 0.70 : substantia nigra (presynaptic) vulnerability
VULN_MSN : 0.75 : striatal medium spiny neuron (POSTsynaptic) vulnerability
VULN_OPC : 0.35 : olivopontocerebellar vulnerability
VULN_ONUF: 1.05 : Onuf nucleus vulnerability
VULN_RESP: 0.40 : arcuate/raphe/pre-Botzinger respiratory vulnerability
VULN_PG : 0.08 : POSTGANGLIONIC vulnerability — near-zero in MSA, high in PAF
COQ2F : 1.0 : CoQ10 biosynthetic capacity (COQ2 variant carriers < 1)
// ---------- glial alpha-synuclein pathology (per-year rates) -------------
KSYN_A : 1.00 : oligodendroglial alpha-syn monomer supply (per yr)
KUPT : 2.20 : uptake of extracellular seed into oligodendrocyte (per yr)
KNUC : 0.85 : p25alpha-templated nucleation (per yr)
KOXA : 0.45 : oxidative promotion of oligomerisation (per yr)
KFIB : 1.10 : oligomer -> fibrillar GCI (per yr)
KCLR : 0.05 : GCI clearance (per yr)
KDEGM : 1.60 : monomer degradation via autophagy/UPS (per yr)
KDEGO : 0.55 : oligomer degradation via autophagy (per yr)
KAUT : 1.30 : oligomer clogging of autophagic flux
KREL : 0.30 : seed release from GCI (per yr)
KRELO : 0.55 : seed release of oligomers (per yr)
KCLRS : 3.20 : glymphatic/microglial seed clearance (per yr)
KP25 : 1.40 : p25alpha relocalisation rate (per yr)
KP25OFF : 0.35 : p25alpha return to myelin (per yr)
HP25 : 0.45 : baseline p25alpha availability floor for nucleation
// ---------- myelin, inflammation, mitochondria ---------------------------
KMYEREP : 0.30 : myelin repair (per yr)
KMYELOSS: 0.30 : myelin loss per unit GCI + oxidant (per yr)
KHOCL : 0.55 : myeloperoxidase/HOCl contribution to myelin loss
KMGL : 2.50 : microglial activation gain (per yr)
KDEB : 0.60 : myelin-debris contribution to microglial activation
KMGLOFF : 2.20 : microglial deactivation (per yr)
KMPO : 1.60 : myeloperoxidase induction by activated microglia (per yr)
KMPOOFF : 1.50 : myeloperoxidase turnover (per yr)
KOX : 1.30 : oxidant generation from MPO (per yr)
KOXCQ : 0.90 : oxidant generation from CoQ10 deficit
KOXOFF : 2.00 : antioxidant handling (per yr)
TAU_CQ : 0.50 : CoQ10 pool time constant (yr)
KCQOX : 0.35 : oxidative consumption of CoQ10
WOXS : 0.55 : weight of oxidative stress in neurodegeneration drive
WINF : 0.40 : weight of neuroinflammation in the drive
WTROPH : 0.85 : weight of LOST OLIGODENDROGLIAL TROPHIC SUPPORT in the drive
KNFL : 55.0 : neurofilament-light release per unit neuronal loss rate
KNFLEL : 0.35 : NfL elimination (per h)
NFL0 : 12.0 : baseline plasma NfL (pg/mL)
// ---------- central autonomic gain --------------------------------------
// The reflex arc is a SERIAL CHAIN (NTS -> CVLM -> RVLM -> IML), so a partial
// loss at each relay multiplies. Exponents > 1 encode that compounding; with
// exponents of 1 a patient who has lost 40% of each relay still retains ~35%
// of reflex gain, which is far too much function to match the clinic.
GEXP_MED: 2.20 : exponent on medullary survival in central gain
GEXP_IML: 1.60 : exponent on IML survival in central gain
FTON : 0.22 : lesion-independent (spinal/local) fraction of tonic drive
SNA0 : 1.00 : reference tonic sympathetic activity
SNAMAX : 2.50 : ceiling on preganglionic sympathetic drive
GBARO : 0.160: baroreflex gain (per mmHg error)
MAP_SET : 92 : baroreflex operating point (mmHg)
TAU_SNA : 0.003: sympathetic activity time constant (h ~ 11 s)
// ---------- postganglionic NE handling ----------------------------------
KREL_NE : 43.0 : NE release rate constant
KA2 : 6.00 : alpha2 autoreceptor affinity (NES units)
KUP : 30.0 : NET reuptake clearance (per h)
KNEDEG : 6.00 : MAO/COMT clearance of synaptic NE (per h)
KVES : 2.40 : vesicular NE repletion (per h)
KUSE : 0.022: vesicular depletion per unit release
KRECYC : 0.035: fraction of reuptaken NE recycled into vesicles
NEPL0 : 250 : supine plasma NE reference (pg/mL) at reference release
KNEREF : 38.5 : healthy reference NE release rate (spillover denominator)
KSPILL : 1.00 : extra plasma spillover per unit NET blockade
KAADC : 0.007: peripheral AADC conversion of droxidopa (per ng/mL)
IC50CBD : 400 : carbidopa concentration for 50% AADC block (ng/mL)
ENET : 0.92 : maximal fractional NET blockade by atomoxetine
EC50ATX : 20 : atomoxetine EC50 for NET blockade (ng/mL)
EYOH : 0.70 : maximal alpha2 blockade by yohimbine
EC50YOH : 55 : yohimbine EC50 (ng/mL)
// ---------- vascular effector ------------------------------------------
KSUPER : 0.60 : denervation supersensitivity gain on alpha1 density
TAU_A1 : 336 : alpha1 receptor up-regulation time constant (h, ~2 wk)
EC50A1 : 1.00 : agonist occupancy EC50 (NES units)
KDMM_EQ : 45 : ng/mL desglymidodrine equivalent to 1 NES unit
ETPR : 1.10 : maximal fractional rise in TPR from alpha1 occupancy
TPR0 : 1.00 : reference peripheral resistance
TAU_TPR : 0.05 : vascular tone time constant (h ~ 3 min)
EWB : 0.28 : water-bolus osmopressor maximal TPR effect
EC50WB : 210 : water-bolus effect-site EC50 (mL equivalent)
KVEN : 2.60 : venous capacitance reduction per unit alpha1 occupancy
KBIND : 0.55 : venous capacitance reduction from abdominal binder
KOCTV : 0.60 : splanchnic constriction from octreotide (max fraction)
EC50OCT : 0.35 : octreotide effect-site EC50 (ug equivalent)
VPOOL : 1.05 : maximal orthostatic splanchnic/dependent pooling (L)
VMEAL : 0.30 : post-prandial additional splanchnic pooling (L)
TAU_POOL: 0.008: pooling time constant (h ~ 30 s)
// ---------- cardiac and haemodynamics ----------------------------------
BV0 : 5.00 : reference blood volume (L)
KBV : 0.55 : blood volume sensitivity to ECF
ECF0 : 14.0 : reference extracellular fluid volume (L)
FCENT : 0.28 : fraction of blood volume that is cardiopulmonary/central
VC0 : 1.40 : reference central blood volume (L) = FCENT*BV0
TAU_VC : 0.006: central volume equilibration (h)
SV0 : 72 : reference stroke volume (mL)
HSV : 0.85 : Frank-Starling exponent on central volume
KINO : 0.12 : sympathetic inotropy gain
HR0 : 66 : intrinsic heart rate reference (bpm)
KCHR : 0.30 : sympathetic chronotropy gain
KVAG : 0.28 : vagal bradycardic gain
VAG0 : 1.00 : reference vagal tone
KVAGP : 0.85 : vagal withdrawal on standing
TAU_VAG : 0.004: vagal time constant (h)
TAU_HR : 0.005: heart-rate time constant (h)
KMAP : 11.1 : MAP per unit cardiac output x resistance (mmHg / (L/min))
TAU_MAP : 0.002: arterial pressure time constant (h ~ 7 s)
KPP : 0.62 : pulse pressure per mL stroke volume (mmHg/mL)
CBFFLOOR: 58 : MAP below which cerebral autoregulation fails (mmHg)
// ---------- renal / volume / endocrine ---------------------------------
NAIN : 6.25 : dietary sodium intake (mmol/h; 150 mmol/day)
WIN : 0.094: water intake (L/h; ~2.25 L/day)
UNA0 : 6.25 : reference sodium excretion (mmol/h)
MAPNAT : 92 : MAP at which natriuresis equals intake (mmHg)
KPN : 0.30 : PRESSURE-natriuresis steepness (per 10 mmHg) — deliberately
KVOLN : 6.0 : VOLUME-error natriuresis steepness (tight homeostasis)
ECFTGT : 14.0 : renal extracellular-volume target (L)
KMR : 0.80 : exponent of mineralocorticoid sodium retention
UV0 : 0.130: reference urine output scale (L/h)
KPNW : 0.22 : pressure-diuresis steepness (per 10 mmHg)
KVOLW : 4.0 : volume-error diuresis steepness
KTHIRST : 1.50 : osmotic thirst gain (fractional water intake per 10% Na rise)
FAQ : 0.62 : maximal fractional water reabsorption from AQP2
KAQ : 0.55 : AQP2 EC50 for water reabsorption
EAQ : 1.00 : maximal AQP2 trafficking
EC50AVP : 2.60 : AVP EC50 for V2R/AQP2 (pg/mL)
KDDA : 3.10 : desmopressin AVP-equivalents per ng/mL
TAU_AQP : 1.20 : AQP2 trafficking time constant (h)
AVP0 : 2.10 : basal plasma AVP (pg/mL)
KOSM : 60.0 : osmotic AVP gain (per fractional plasma-sodium change)
KBAVP : 0.16 : BARORECEPTOR AVP gain (per mmHg) — gated by G_CENT
TAU_AVP : 0.35 : AVP time constant (h)
PRA0 : 1.00 : reference plasma renin activity
FRTON : 0.30 : sympathetic-independent renin fraction
KBRENIN : 0.70 : beta1-sympathetic renin gain
KPRAP : 1.30 : renal perfusion-pressure renin gain
TAU_PRA : 0.60 : renin time constant (h)
ALD0 : 1.00 : reference aldosterone
TAU_ALD : 1.50 : aldosterone time constant (h)
EFLU : 2.30 : maximal fludrocortisone mineralocorticoid effect
EC50FLU : 1.40 : fludrocortisone effect-site EC50 (ng/mL)
EPYR : 0.85 : maximal pyridostigmine amplification of BAROREFLEX traffic
EC50PYR : 26 : pyridostigmine EC50 (ng/mL)
// ---------- dopaminergic ------------------------------------------------
KDAEND : 1.00 : endogenous striatal dopamine per unit nigral survival
KDALD : 0.85 : dopamine produced per ug/mL brain levodopa
FGLIA : 0.30 : fraction of levodopa decarboxylation independent of terminals
EC50DA : 0.55 : striatal dopamine EC50 for postsynaptic signalling
GEXP_MSN: 1.15 : exponent on MSN survival in postsynaptic gain
TAU_DA : 0.30 : striatal dopamine time constant (h)
PMAX : 26 : maximal parkinsonism contribution to UMSARS-II
AMAX : 22 : maximal ataxia contribution to UMSARS-II
PVRMAX : 320 : maximal post-void residual (mL)
SMAX : 1.00 : maximal stridor severity index
// ---------- clinical endpoints -----------------------------------------
SBPTH : 96 : standing SBP at which OH symptoms are half-maximal (mmHg)
OHSLOPE : 9.0 : steepness of the symptom-vs-standing-SBP curve (mmHg)
OHMAX : 10.0 : maximal OHSA item-1 style score
TAU_OHSA: 24 : symptom-score time constant (h)
TAU_U : 720 : UMSARS time constant (h, ~1 month)
WPARK : 1.00 : weight of parkinsonism in UMSARS-II
WATAX : 1.00 : weight of ataxia in UMSARS-II
WBULB : 6.0 : weight of bulbar/respiratory involvement in UMSARS-II
WOH : 0.55 : weight of OH symptoms in UMSARS-I
WURO : 3.5 : weight of urogenital failure in UMSARS-I
WS_IML : 12 : SCOPA-AUT weight, sympathetic failure
WS_URO : 10 : SCOPA-AUT weight, urogenital
WS_GI : 8 : SCOPA-AUT weight, gastrointestinal (via DMV proxy NRESP)
WS_OH : 1.4 : SCOPA-AUT weight, orthostatic symptoms
H0 : 0.0000017 : baseline hazard (per h)
B_U2 : 0.42 : log-hazard per 10 UMSARS-II points
B_STR : 0.85 : log-hazard for stridor
B_OH : 0.22 : log-hazard per 5 OHSA points
B_RESP : 1.15 : log-hazard for respiratory nucleus loss
// ---------- drug pharmacokinetics --------------------------------------
KA_MID : 3.00 : midodrine absorption (per h)
KM_MID : 1.85 : midodrine -> desglymidodrine conversion (per h)
KE_MID : 0.45 : midodrine non-converting elimination (per h)
VMID : 60 : midodrine central volume (L)
KE_DMM : 0.231: desglymidodrine elimination (per h; t1/2 3 h)
VDMM : 95 : desglymidodrine volume (L)
KA_DRX : 0.90 : droxidopa absorption (per h)
KE_DRX : 0.277: droxidopa elimination (per h; t1/2 2.5 h)
VDRX : 35 : droxidopa volume (L)
KE_FLU : 0.029: fludrocortisone effective elimination (per h; t1/2 24 h)
VFLU : 45 : fludrocortisone volume (L)
KA_FLU : 1.60 : fludrocortisone absorption (per h)
KE_ATX : 0.133: atomoxetine elimination (per h; t1/2 5.2 h extensive metaboliser)
CYP2D6 : 1.0 : CYP2D6 activity multiplier (0.24 = poor metaboliser)
VATX : 62 : atomoxetine volume (L)
KA_ATX : 2.00 : atomoxetine absorption (per h)
KE_PYR : 0.198: pyridostigmine elimination (per h; t1/2 3.5 h)
VPYR : 70 : pyridostigmine apparent volume (L)
KA_PYR : 1.00 : pyridostigmine absorption (per h)
KE_DDA : 0.231: desmopressin elimination (per h; t1/2 3 h)
VDDA : 25 : desmopressin volume (L)
KE_LD : 0.462: levodopa elimination (per h; t1/2 1.5 h with carbidopa)
VLD : 40 : levodopa central volume (L)
KA_LD : 2.00 : levodopa absorption (per h, gastroparesis-gated)
KGP : 0.55 : fractional reduction of KA_LD by gastroparesis
KIN_LDB : 0.35 : levodopa brain influx (LAT1) (per h)
KOUT_LDB: 0.95 : brain levodopa efflux + decarboxylation (per h)
VLDB : 1.40 : brain levodopa distribution volume (L)
KE_CBD : 0.347: carbidopa elimination (per h; t1/2 2 h)
VCBD : 50 : carbidopa volume (L)
KA_CBD : 1.50 : carbidopa absorption (per h)
KE_UBQ : 0.020: ubiquinol elimination (per h; t1/2 ~34 h)
VUBQ : 120 : ubiquinol volume (L)
KA_UBQ : 0.25 : ubiquinol absorption (per h, low bioavailability)
EUBQ : 0.55 : maximal CoQ10 repletion effect
EC50UBQ : 2200 : ubiquinol EC50 (ng/mL)
KE_VRD : 0.058: verdiperstat elimination (per h; t1/2 12 h)
VVRD : 55 : verdiperstat volume (L)
KA_VRD : 1.20 : verdiperstat absorption (per h)
EVRD : 0.85 : maximal myeloperoxidase inhibition
EC50VRD : 900 : verdiperstat EC50 (ng/mL)
KE_MAB : 0.00138 : anti-alpha-syn mAb elimination (per h; t1/2 ~21 d)
VMAB : 3.60 : antibody central volume (L)
KMABSEED: 0.42 : seed neutralisation per mg/L antibody (per yr)
KE_YOH : 0.462: yohimbine elimination (per h)
VYOH : 80 : yohimbine volume (L)
KA_YOH : 2.20 : yohimbine absorption (per h)
KE_OCT : 0.385: octreotide elimination (per h; t1/2 1.8 h)
VOCT : 20 : octreotide volume (L)
KE_WB : 0.92 : water-bolus effect-site decay (per h)
IASO : 0.0 : SNCA-lowering ASO fractional knockdown (0-0.8)
GLYMPH : 0.0 : additional glymphatic seed clearance (0-1)
// ---------- multiscale housekeeping ------------------------------------
// TSCALE stretches ONLY the fast cardiovascular time constants (seconds to
// minutes). The disease layer, the renal/endocrine layer and every
// quasi-steady VALUE are untouched, so a multi-year run with TSCALE = 40
// (TAU_MAP 7 s -> 5 min) reproduces the same trajectory as TSCALE = 1 at a
// small fraction of the cost: a 61-state stiff system with a numerical
// Jacobian otherwise spends essentially all of its time resolving a
// 7-second transient that no annual endpoint depends on. Tilt runs and
// 48-h runs use TSCALE = 1, where those seconds ARE the measurement.
// msa_selftest() asserts the two agree to better than 1%.
TSCALE : 1 : stretch factor on the fast haemodynamic time constants
// INITMAIN gates the explicit compartment initialisation in $MAIN. It must
// be 1 for a fresh healthy start and 0 whenever a run is RESTARTED from a
// saved state via init(): mrgsolve applies $MAIN after init(), so leaving
// the X_0 assignments unguarded silently resets a diseased patient back to
// healthy physiology and every tilt/48-h experiment then reports the SAME
// numbers regardless of disease stage or phenotype.
INITMAIN: 1 : 1 = initialise from parameters, 0 = keep init() state
HPY : 8766 : hours per year (converts per-year rates to per-hour)
$CMT @annotated
// --- glial proteinopathy (10) ---
P25A : p25alpha/TPPP mislocalised to soma (fraction)
ASYNM : oligodendroglial alpha-syn monomer (a.u.)
ASYNO : alpha-syn oligomer (a.u.)
GCI : glial cytoplasmic inclusion burden (a.u.)
SEED : extracellular seeding-competent alpha-syn (a.u.)
MGL : activated microglia (1 = resting baseline)
MPO : myeloperoxidase activity (a.u.)
OXS : oxidative stress (a.u.)
CQ : coenzyme Q10 pool (1 = normal)
MYE : myelin integrity (1 = normal)
// --- neuronal populations, surviving fraction (9) ---
NIML : spinal IML preganglionic sympathetic neurons
NMED : NTS/CVLM/RVLM medullary cardiovascular neurons
NVAG : cardiovagal preganglionic neurons
NSN : substantia nigra dopaminergic neurons (PREsynaptic)
NMSN : striatal medium spiny neurons (POSTsynaptic)
NOPC : olivopontocerebellar neurons
NONUF : Onuf nucleus motor neurons
NRESP : arcuate/raphe/pre-Botzinger respiratory neurons
POSTG : POSTGANGLIONIC sympathetic terminal integrity
// --- biomarker (1) ---
NFL : plasma neurofilament light (pg/mL)
// --- cardiovascular fast loop (10) ---
SNA : sympathetic nerve activity (1 = healthy supine)
NES : synaptic norepinephrine at vascular alpha1 (a.u.)
NEVES : releasable vesicular NE pool (1 = full)
A1R : vascular alpha1-adrenoceptor density (1 = normal)
TPR : total peripheral resistance (1 = reference)
VSP : splanchnic/dependent pooled volume (L)
VC : central blood volume (L)
HR : heart rate (bpm)
VAG : cardiovagal tone (1 = reference)
MAP : mean arterial pressure (mmHg)
// --- renal / volume / endocrine (6) ---
NAB : total body exchangeable sodium (mmol)
ECF : extracellular fluid volume (L)
AQP2 : collecting-duct AQP2 apical trafficking (1 = maximal)
AVP : plasma arginine vasopressin (pg/mL)
PRA : plasma renin activity (1 = reference)
ALD : aldosterone (1 = reference)
// --- dopaminergic (1) ---
DA : striatal synaptic dopamine (a.u.)
// --- drug PK (19) ---
MIDG : midodrine gut (mg)
MIDC : midodrine plasma (mg)
DMM : desglymidodrine plasma (mg)
DRXG : droxidopa gut (mg)
DRXC : droxidopa plasma (mg)
FLUG : fludrocortisone gut (mg)
FLU : fludrocortisone plasma (mg)
ATXG : atomoxetine gut (mg)
ATX : atomoxetine plasma (mg)
PYRG : pyridostigmine gut (mg)
PYR : pyridostigmine plasma (mg)
DDA : desmopressin plasma (ug)
LDG : levodopa gut (mg)
LDC : levodopa plasma (mg)
LDB : brain levodopa (mg)
CBDG : carbidopa gut (mg)
CBD : carbidopa plasma (mg)
UBQG : ubiquinol gut (mg)
UBQ : ubiquinol plasma (mg)
// --- drug PK, continued (5) ---
VRDG : verdiperstat gut (mg)
VRD : verdiperstat plasma (mg)
MAB : anti-alpha-syn antibody plasma (mg)
YOHG : yohimbine gut (mg)
YOH : yohimbine plasma (mg)
OCT : octreotide plasma (ug)
WB : water-bolus effect site (mL equivalent)
// --- clinical endpoints (5) ---
OHSA : orthostatic-hypotension symptom score (0-10)
U1 : UMSARS Part I (activities of daily living)
U2 : UMSARS Part II (motor examination)
SCOPA : SCOPA-AUT autonomic symptom score
CUMH : cumulative mortality hazard
$MAIN
// ---- initial conditions: healthy physiology, zero pathology ----
// Guarded by INITMAIN so that init() can restart from a frozen disease state.
if(INITMAIN > 0.5) {
P25A_0 = HP25;
ASYNM_0 = KSYN_A/KDEGM;
ASYNO_0 = 0.0;
GCI_0 = 0.0;
SEED_0 = 0.0;
MGL_0 = 1.0;
MPO_0 = 1.0;
OXS_0 = 1.0;
CQ_0 = COQ2F;
MYE_0 = 1.0;
NIML_0 = 1.0; NMED_0 = 1.0; NVAG_0 = 1.0; NSN_0 = 1.0;
NMSN_0 = 1.0; NOPC_0 = 1.0; NONUF_0= 1.0; NRESP_0 = 1.0;
POSTG_0= 1.0;
NFL_0 = NFL0;
SNA_0 = SNA0;
NES_0 = 1.0;
NEVES_0 = 1.0;
A1R_0 = 1.0;
TPR_0 = TPR0;
VSP_0 = 0.0;
VC_0 = FCENT*BV0;
HR_0 = HR0;
VAG_0 = VAG0;
MAP_0 = MAP_SET;
ECF_0 = ECF0;
NAB_0 = 140.0*ECF0;
AQP2_0 = 0.45;
AVP_0 = AVP0;
PRA_0 = PRA0;
ALD_0 = ALD0;
DA_0 = KDAEND;
OHSA_0 = 0.0;
U1_0 = 0.0;
U2_0 = 0.0;
SCOPA_0 = 0.0;
CUMH_0 = 0.0;
}
F_LDG = 1.0;
ALAG_LDG = 0.0;
$ODE
// =====================================================================
// 0. posture / meal drivers (smooth, so the solver never sees a step)
// =====================================================================
double tod = fmod(SOLVERTIME, 24.0);
double wake = SIGM(tod - WAKE_H, 8.0) * SIGM(SLEEP_H - tod, 8.0);
double POS;
if(POSMODE < 0.5) POS = 0.0;
else if(POSMODE < 1.5) POS = 1.0;
else if(POSMODE < 2.5) POS = UPFRAC*wake + HUTF*(1.0 - wake);
else POS = SIGM(SOLVERTIME - TILT_T, 120.0);
POS = POS*(1.0 - 0.25*COUNTERM); // counter-manoeuvres offload the pool
double meal = MEALAMP*wake*MEALB(tod);
// =====================================================================
// 1. glial alpha-synuclein pathology (rates given per year -> per hour)
// =====================================================================
double AUTOF = 1.0/(1.0 + KAUT*ASYNO); // autophagic flux, clogged by oligomers
double p25 = HP25 + (1.0 - HP25)*P25A; // effective nucleation template
double CMAB = MAB/VMAB; // mg/L
dxdt_P25A = (KP25*(1.0 - MYE)*(1.0 - P25A) - KP25OFF*P25A)/HPY;
dxdt_ASYNM = (KSYN_A*(1.0 - IASO) + KUPT*SEED
- KNUC*ASYNM*p25 - KOXA*OXS*ASYNM*0.5
- KDEGM*ASYNM*AUTOF)/HPY;
dxdt_ASYNO = (KNUC*ASYNM*p25 + KOXA*OXS*ASYNM*0.5
- KFIB*ASYNO - KDEGO*ASYNO*AUTOF)/HPY;
dxdt_GCI = (KFIB*ASYNO - KCLR*GCI)/HPY;
dxdt_SEED = (KREL*GCI + KRELO*ASYNO
- KCLRS*(1.0 + GLYMPH)*SEED
- KMABSEED*CMAB*SEED)/HPY;
// =====================================================================
// 2. myelin, neuroinflammation, mitochondria
// =====================================================================
double MGLX = (MGL > 1.0) ? (MGL - 1.0) : 0.0;
double CUBQ = 1000.0*UBQ/VUBQ; // ng/mL
double CVRD = 1000.0*VRD/VVRD; // ng/mL
double IVRD = EVRD*CVRD/(EC50VRD + CVRD);
double CQss = COQ2F*(1.0 + EUBQ*CUBQ/(EC50UBQ + CUBQ))/(1.0 + KCQOX*(OXS - 1.0 > 0 ? OXS - 1.0 : 0.0));
dxdt_MGL = (KMGL*(SEED + KDEB*(1.0 - MYE)) - KMGLOFF*MGLX)/HPY;
dxdt_MPO = (KMPO*MGLX*(1.0 - IVRD) - KMPOOFF*(MPO - 1.0))/HPY;
dxdt_OXS = (KOX*(MPO - 1.0) + KOXCQ*(1.0 - CQ) - KOXOFF*(OXS - 1.0))/HPY;
dxdt_CQ = (CQss - CQ)/(TAU_CQ*HPY);
dxdt_MYE = (KMYEREP*(1.0 - MYE) - KMYELOSS*(GCI + KHOCL*(MPO - 1.0))*MYE)/HPY;
// =====================================================================
// 3. neurodegeneration — ONE drive, many vulnerabilities
// The trophic term is what makes MSA a GLIAL disease: neurons die
// because their oligodendrocytes stopped supporting them.
// =====================================================================
double NDRIVE = KND*(GCI + WOXS*(OXS - 1.0) + WINF*MGLX + WTROPH*(1.0 - MYE))/HPY;
if(NDRIVE < 0.0) NDRIVE = 0.0;
dxdt_NIML = -NDRIVE*VULN_IML *NIML;
dxdt_NMED = -NDRIVE*VULN_MED *NMED;
dxdt_NVAG = -NDRIVE*VULN_VAG *NVAG;
dxdt_NSN = -NDRIVE*VULN_SN *NSN;
dxdt_NMSN = -NDRIVE*VULN_MSN *NMSN;
dxdt_NOPC = -NDRIVE*VULN_OPC *NOPC;
dxdt_NONUF = -NDRIVE*VULN_ONUF*NONUF;
dxdt_NRESP = -NDRIVE*VULN_RESP*NRESP;
dxdt_POSTG = -NDRIVE*VULN_PG *POSTG;
double LOSSRATE = NDRIVE*(VULN_IML*NIML + VULN_MED*NMED + VULN_SN*NSN
+ VULN_MSN*NMSN + VULN_OPC*NOPC + VULN_RESP*NRESP)/6.0;
dxdt_NFL = KNFL*LOSSRATE*HPY - KNFLEL*(NFL - NFL0);
// =====================================================================
// 4. THE PRODUCT — central gain vs postganglionic integrity
// =====================================================================
double G_CENT = pow(NMED > 1e-6 ? NMED : 1e-6, GEXP_MED)
* pow(NIML > 1e-6 ? NIML : 1e-6, GEXP_IML);
double CPYR = 1000.0*PYR/VPYR;
double PYRAMP= 1.0 + EPYR*CPYR/(EC50PYR + CPYR); // amplifies REFLEX traffic only
double berr = MAP_SET - MAP;
double SNAss = SNA0*(FTON + (1.0 - FTON)*G_CENT) + GBARO*G_CENT*PYRAMP*berr;
if(SNAss < 0.02) SNAss = 0.02;
if(SNAss > SNAMAX) SNAss = SNAMAX; // preganglionic drive cannot rise without limit
dxdt_SNA = (SNAss - SNA)/(TAU_SNA*TSCALE);
// =====================================================================
// 5. norepinephrine handling at the postganglionic terminal
// =====================================================================
double CYOH = 1000.0*YOH/VYOH;
double IYOH = EYOH*CYOH/(EC50YOH + CYOH);
double A2INH = 1.0/(1.0 + (NES/KA2)*(1.0 - IYOH));
double NEREL = KREL_NE*SNA*NEVES*POSTG*A2INH;
double CDRX = 1000.0*DRXC/VDRX;
double CCBD = 1000.0*CBD/VCBD;
double CARBI = CCBD/(IC50CBD + CCBD); // <-- shared AADC term
double DROXNE= KAADC*CDRX*POSTG*(1.0 - CARBI); // needs terminals AND no carbidopa
double CATX = 1000.0*ATX/VATX;
double NETI = ENET*CATX/(EC50ATX + CATX);
// NET is expressed ON the postganglionic terminal, so reuptake capacity is
// gated by POSTG. This is what makes NET inhibition phenotype-selective:
// where the terminal is gone there is no transporter left to block and
// atomoxetine has nothing to do. Without the gate a PAF patient gets almost
// the same atomoxetine response as an MSA patient, which is the opposite of
// what the crossover trials show.
double NEUP = KUP*POSTG*(1.0 - NETI)*NES;
dxdt_NES = NEREL + DROXNE - NEUP - KNEDEG*NES;
dxdt_NEVES = KVES*POSTG*(1.0 - NEVES) - KUSE*NEREL + KRECYC*NEUP;
// =====================================================================
// 6. vascular effector — denervation supersensitivity multiplies EVERY
// alpha1 agonist, endogenous or pharmacological
// =====================================================================
double A1ss = 1.0 + KSUPER*(1.0 - G_CENT);
dxdt_A1R = (A1ss - A1R)/TAU_A1;
double CDMM = 1000.0*DMM/VDMM;
double AGON_N = NES;
double AGON = NES + CDMM/KDMM_EQ;
double OCC = AGON/(EC50A1 + AGON);
double WBEF = EWB*WB/(EC50WB + WB);
double TPRss= TPR0*(1.0 + ETPR*A1R*OCC + WBEF);
dxdt_TPR = (TPRss - TPR)/(TAU_TPR*TSCALE);
double COCT = OCT/VOCT;
double OCTEF = KOCTV*COCT/(EC50OCT + COCT);
double CAPREL = 1.0/(1.0 + KVEN*A1R*OCC + KBIND*BINDER);
double POOLss = (POS*VPOOL + meal*VMEAL*(1.0 - OCTEF))*CAPREL*(1.0 - 0.55*OCTEF*0.0);
dxdt_VSP = (POOLss - VSP)/(TAU_POOL*TSCALE);
// =====================================================================
// 7. cardiac output and arterial pressure
// =====================================================================
double BV = BV0*(1.0 + KBV*(ECF/ECF0 - 1.0));
double VCss = FCENT*BV - VSP;
if(VCss < 0.15) VCss = 0.15;
dxdt_VC = (VCss - VC)/(TAU_VC*TSCALE);
double SV_ = SV0*pow(VC/VC0 > 0.05 ? VC/VC0 : 0.05, HSV)*(1.0 + KINO*SNA*POSTG);
double GMED = pow(NMED > 1e-6 ? NMED : 1e-6, GEXP_MED);
double VAGss= VAG0*pow(NVAG > 1e-6 ? NVAG : 1e-6, 1.0)*(1.0 - KVAGP*POS*GMED);
dxdt_VAG = (VAGss - VAG)/(TAU_VAG*TSCALE);
double HRss = HR0*(1.0 + KCHR*SNA*POSTG)/(1.0 + KVAG*VAG);
dxdt_HR = (HRss - HR)/(TAU_HR*TSCALE);
double CO_ = SV_*HR/1000.0; // L/min
double MAPss= KMAP*CO_*TPR;
dxdt_MAP = (MAPss - MAP)/(TAU_MAP*TSCALE);
// =====================================================================
// 8. renal sodium / water — THE NOCTURNAL LOOP
// Both natriuresis and diuresis rise exponentially with MAP, so supine
// hypertension at night spends the very volume the patient needs at 07:00.
// =====================================================================
double CFLU = 1000.0*FLU/VFLU;
double MR = ALD/ALD0 + EFLU*CFLU/(EC50FLU + CFLU);
if(MR < 0.05) MR = 0.05;
// Sodium excretion carries TWO terms with deliberately different gains:
// a STEEP volume-error term (long-run homeostasis, keeps total body sodium
// bounded) and a SHALLOW pressure term (~1.5-2x over a 25 mmHg supine rise).
// It is the shallow term that closes the nocturnal loop; making it steep
// would make the loop implausibly violent, making it zero would abolish the
// central clinical trade-off of the disease.
double VOLERR = ECF/ECFTGT;
if(VOLERR < 0.3) VOLERR = 0.3;
double UNAV = UNA0*pow(VOLERR, KVOLN)*exp(KPN*(MAP - MAPNAT)/10.0)/pow(MR, KMR);
dxdt_NAB = NAIN - UNAV;
if(NAB < 200.0 && dxdt_NAB < 0.0) dxdt_NAB = 0.0;
double CDDA = 1000.0*DDA/VDDA; // pg/mL-equivalent scale
double AQPss= EAQ*(AVP + KDDA*CDDA)/(EC50AVP + AVP + KDDA*CDDA);
dxdt_AQP2 = (AQPss - AQP2)/TAU_AQP;
double UVOL = UV0*pow(VOLERR, KVOLW)*exp(KPNW*(MAP - MAPNAT)/10.0)
*(1.0 - FAQ*AQP2/(KAQ + AQP2));
if(UVOL < 0.005) UVOL = 0.005;
double NACONC = NAB/(ECF > 1.0 ? ECF : 1.0);
// osmotic thirst — the intact half of volume defence in MSA
double WINE = WIN*(1.0 + KTHIRST*(NACONC/140.0 - 1.0)*10.0);
if(WINE < 0.3*WIN) WINE = 0.3*WIN;
dxdt_ECF = WINE - UVOL;
double AVPss = AVP0*(1.0 + KOSM*(NACONC/140.0 - 1.0))
+ KBAVP*G_CENT*(berr > 0.0 ? berr : 0.0); // baroreflex arm is GATED
if(AVPss > 40.0) AVPss = 40.0;
if(AVPss < 0.1) AVPss = 0.1;
dxdt_AVP = (AVPss - AVP)/TAU_AVP;
double PRAss = PRA0*(FRTON + KBRENIN*SNA*POSTG)
*pow(MAPNAT/(MAP > 20.0 ? MAP : 20.0), KPRAP);
dxdt_PRA = (PRAss - PRA)/TAU_PRA;
double ALDss = ALD0*pow(PRA/PRA0 > 0.05 ? PRA/PRA0 : 0.05, 0.8);
dxdt_ALD = (ALDss - ALD)/TAU_ALD;
// =====================================================================
// 9. dopaminergic: PREsynaptic supply vs POSTsynaptic target
// =====================================================================
double CLDB = LDB/VLDB; // mg/L = ug/mL
double DAss = KDAEND*NSN + KDALD*CLDB*(FGLIA + (1.0 - FGLIA)*NSN);
dxdt_DA = (DAss - DA)/TAU_DA;
// =====================================================================
// 10. drug pharmacokinetics
// =====================================================================
dxdt_MIDG = -KA_MID*MIDG;
dxdt_MIDC = KA_MID*MIDG - (KM_MID + KE_MID)*MIDC;
dxdt_DMM = KM_MID*MIDC - KE_DMM*DMM;
dxdt_DRXG = -KA_DRX*DRXG;
dxdt_DRXC = KA_DRX*DRXG - KE_DRX*DRXC;
dxdt_FLUG = -KA_FLU*FLUG;
dxdt_FLU = KA_FLU*FLUG - KE_FLU*FLU;
dxdt_ATXG = -KA_ATX*ATXG;
dxdt_ATX = KA_ATX*ATXG - (KE_ATX*CYP2D6)*ATX;
dxdt_PYRG = -KA_PYR*PYRG;
dxdt_PYR = KA_PYR*PYRG - KE_PYR*PYR;
dxdt_DDA = -KE_DDA*DDA;
double KA_LDe = KA_LD*(1.0 - KGP*(1.0 - NRESP)); // gastroparesis gates absorption
dxdt_LDG = -KA_LDe*LDG;
dxdt_LDC = KA_LDe*LDG - KE_LD*LDC - KIN_LDB*LDC;
dxdt_LDB = KIN_LDB*LDC - KOUT_LDB*LDB;
dxdt_CBDG = -KA_CBD*CBDG;
dxdt_CBD = KA_CBD*CBDG - KE_CBD*CBD;
dxdt_UBQG = -KA_UBQ*UBQG;
dxdt_UBQ = KA_UBQ*UBQG - KE_UBQ*UBQ;
dxdt_VRDG = -KA_VRD*VRDG;
dxdt_VRD = KA_VRD*VRDG - KE_VRD*VRD;
dxdt_MAB = -KE_MAB*MAB;
dxdt_YOHG = -KA_YOH*YOHG;
dxdt_YOH = KA_YOH*YOHG - KE_YOH*YOH;
dxdt_OCT = -KE_OCT*OCT;
dxdt_WB = -KE_WB*WB;
// =====================================================================
// 11. clinical endpoints
// =====================================================================
// Plasma norepinephrine is measured SPILLOVER — a function of RELEASE and of
// how much escapes reuptake — not of the synaptic concentration that drives
// the receptor. Reading it off NES makes a phenotype with dead terminals but
// intact central drive look near-normal, inverting the single most useful
// bedside discriminator between MSA (normal supine NE, blunted increment) and
// PAF (low supine NE).
double NEPLA = NEPL0*(NEREL/KNEREF)*(1.0 + KSPILL*NETI);
double PPr = KPP*SV_;
double SBPn = MAP + PPr/2.0;
double OHSt = OHMAX/(1.0 + exp((SBPn - SBPTH)/OHSLOPE));
dxdt_OHSA = (OHSt - OHSA)/TAU_OHSA;
double G_POST = pow(NMSN > 1e-6 ? NMSN : 1e-6, GEXP_MSN);
double STRIAT = G_POST*DA/(EC50DA + DA);
double STR0 = KDAEND/(EC50DA + KDAEND);
double PARK = PMAX*(1.0 - STRIAT/STR0);
if(PARK < 0.0) PARK = 0.0;
double ATAX = AMAX*(1.0 - NOPC);
double PVRES = PVRMAX*(1.0 - NONUF);
double STRID = SMAX*(1.0 - NRESP);
double U2t = WPARK*PARK + WATAX*ATAX + WBULB*STRID;
dxdt_U2 = (U2t - U2)/TAU_U;
double U1t = 0.80*U2 + WOH*OHSA + WURO*PVRES/PVRMAX*4.0;
dxdt_U1 = (U1t - U1)/TAU_U;
double SCt = WS_IML*(1.0 - NIML) + WS_URO*(1.0 - NONUF)
+ WS_GI*(1.0 - NRESP) + WS_OH*OHSA;
dxdt_SCOPA = (SCt - SCOPA)/TAU_U;
double HAZ = H0*exp(B_U2*U2/10.0 + B_STR*STRID + B_OH*OHSA/5.0 + B_RESP*(1.0 - NRESP));
dxdt_CUMH = HAZ;
$TABLE
// ---- RECOMPUTED at the exact output time --------------------------------
// mrgsolve hoists every `double` declaration from $ODE/$TABLE into one
// shared scope, so the algebra below deliberately RE-ASSIGNS (no `double`)
// the quantities that $ODE also computes. That matters: the last derivative
// evaluation the solver makes is not generally at the output time, and with
// TAU_MAP = 0.002 h a stale MAP-derived value would corrupt the tilt test.
tod = fmod(TIME, 24.0);
wake = SIGM(tod - WAKE_H, 8.0)*SIGM(SLEEP_H - tod, 8.0);
double POSo;
if(POSMODE < 0.5) POSo = 0.0;
else if(POSMODE < 1.5) POSo = 1.0;
else if(POSMODE < 2.5) POSo = UPFRAC*wake + HUTF*(1.0 - wake);
else POSo = SIGM(TIME - TILT_T, 120.0);
POSo = POSo*(1.0 - 0.25*COUNTERM);
G_CENT = pow(NMED > 1e-6 ? NMED : 1e-6, GEXP_MED)
* pow(NIML > 1e-6 ? NIML : 1e-6, GEXP_IML);
G_POST = pow(NMSN > 1e-6 ? NMSN : 1e-6, GEXP_MSN);
SV_ = SV0*pow(VC/VC0 > 0.05 ? VC/VC0 : 0.05, HSV)*(1.0 + KINO*SNA*POSTG);
CO_ = SV_*HR/1000.0;
PPr = KPP*SV_;
double SBP = MAP + PPr/2.0;
double DBP = MAP - PPr/3.0;
CDMM = 1000.0*DMM/VDMM;
CDRX = 1000.0*DRXC/VDRX;
CATX = 1000.0*ATX/VATX;
CCBD = 1000.0*CBD/VCBD;
CFLU = 1000.0*FLU/VFLU;
CPYR = 1000.0*PYR/VPYR;
CLDB = LDB/VLDB;
double CLDC = 1000.0*LDC/VLD;
CUBQ = 1000.0*UBQ/VUBQ;
CVRD = 1000.0*VRD/VVRD;
CMAB = MAB/VMAB;
AGON = NES + CDMM/KDMM_EQ;
OCC = AGON/(EC50A1 + AGON);
CARBI= CCBD/(IC50CBD + CCBD);
NETI = ENET*CATX/(EC50ATX + CATX);
DROXNE = KAADC*CDRX*POSTG*(1.0 - CARBI);
MR = ALD/ALD0 + EFLU*CFLU/(EC50FLU + CFLU);
if(MR < 0.05) MR = 0.05;
VOLERR = ECF/ECFTGT;
if(VOLERR < 0.3) VOLERR = 0.3;
UNAV = UNA0*pow(VOLERR, KVOLN)*exp(KPN*(MAP - MAPNAT)/10.0)/pow(MR, KMR);
UVOL = UV0*pow(VOLERR, KVOLW)*exp(KPNW*(MAP - MAPNAT)/10.0)
*(1.0 - FAQ*AQP2/(KAQ + AQP2));
double NAPL = NAB/(ECF > 1.0 ? ECF : 1.0);
AGON_N = NES;
double IYOHt = EYOH*(1000.0*YOH/VYOH)/(EC50YOH + 1000.0*YOH/VYOH);
double A2INHt = 1.0/(1.0 + (NES/KA2)*(1.0 - IYOHt));
NEREL = KREL_NE*SNA*NEVES*POSTG*A2INHt;
double NEPL = NEPL0*(NEREL/KNEREF)*(1.0 + KSPILL*NETI);
STRIAT = G_POST*DA/(EC50DA + DA);
STR0 = KDAEND/(EC50DA + KDAEND);
PARK = PMAX*(1.0 - STRIAT/STR0);
if(PARK < 0.0) PARK = 0.0;
ATAX = AMAX*(1.0 - NOPC);
PVRES = PVRMAX*(1.0 - NONUF);
STRID = SMAX*(1.0 - NRESP);
double U4 = 1.0 + 4.0*(U2/(U2 + 28.0));
double SURV = exp(-CUMH);
double CBFMARG = MAP - CBFFLOOR;
double YEARS = TIME/8766.0;
$CAPTURE @annotated
YEARS : disease years since pathology onset
POSo : orthostatic load (0 supine, 1 fully upright)
SBP : systolic blood pressure (mmHg)
DBP : diastolic blood pressure (mmHg)
CO_ : cardiac output (L/min)
SV_ : stroke volume (mL)
G_CENT : CENTRAL autonomic gain (product of medullary x IML survival)
G_POST : POSTsynaptic striatal gain
OCC : vascular alpha1 occupancy
NETI : fractional NET blockade by atomoxetine
CARBI : fractional peripheral AADC blockade by carbidopa
DROXNE : NE generated from droxidopa (a.u./h)
NEPL : plasma norepinephrine, spillover-based (pg/mL)
NEREL : postganglionic NE release rate (a.u./h)
MR : mineralocorticoid receptor activation (1 = reference)
UNAV : urinary sodium excretion (mmol/h)
NAPL : plasma sodium concentration (mmol/L)
UVOL : urine output (L/h)
PARK : parkinsonism severity (UMSARS-II points)
ATAX : ataxia severity (UMSARS-II points)
STRIAT : striatal output (postsynaptic gain x dopamine occupancy)
PVRES : post-void residual (mL)
STRID : stridor severity index
U4 : UMSARS Part IV global disability (1-5)
SURV : survival probability
CBFMARG: cerebral autoregulation margin (mmHg above floor)
CDMM : desglymidodrine (ng/mL)
CDRX : droxidopa (ng/mL)
CATX : atomoxetine (ng/mL)
CLDC : plasma levodopa (ng/mL)
CLDB : brain levodopa (ug/mL)
CFLU : fludrocortisone (ng/mL)
CPYR : pyridostigmine (ng/mL)
CUBQ : ubiquinol (ng/mL)
CVRD : verdiperstat (ng/mL)
CMAB : anti-alpha-syn antibody (mg/L)
'
mod <- mcode("msa_qsp", msa_code, soloc = tempdir())
## =====================================================================
## PHENOTYPES — set ONLY by regional vulnerability, never by equations
## =====================================================================
## MSA-P striatonigral predominant (about 60-70% in Western series)
## MSA-C olivopontocerebellar predominant (predominant in Japanese series)
## PAF pure autonomic failure — the postganglionic mirror image
## PD Parkinson disease — presynaptic loss, POSTsynaptic target intact
pheno <- list(
MSA_P = c(VULN_SN = 0.70, VULN_MSN = 0.75, VULN_OPC = 0.35,
VULN_IML = 1.00, VULN_MED = 0.85, VULN_PG = 0.08),
MSA_C = c(VULN_SN = 0.22, VULN_MSN = 0.20, VULN_OPC = 1.05,
VULN_IML = 0.95, VULN_MED = 0.80, VULN_PG = 0.08),
PAF = c(VULN_SN = 0.05, VULN_MSN = 0.05, VULN_OPC = 0.05,
VULN_IML = 0.10, VULN_MED = 0.08, VULN_PG = 1.90,
VULN_ONUF = 0.10, VULN_RESP = 0.05, KND = 0.34),
PD = c(VULN_SN = 1.10, VULN_MSN = 0.02, VULN_OPC = 0.05,
VULN_IML = 0.28, VULN_MED = 0.10, VULN_PG = 0.55,
VULN_ONUF = 0.12, VULN_RESP = 0.10, KND = 0.30)
)
## An empty list is a legal "no override" — mrgsolve rejects it, so gate it.
pp <- function(m, p) if (length(p)) mrgsolve::param(m, p) else m
## =====================================================================
## STAGE 1 — natural history (years). Circadian posture, no drug.
## =====================================================================
msa_history <- function(years = 12, phenotype = "MSA_P", p = list(), delta = 24,
tscale = 40) {
mod %>%
param(as.list(pheno[[phenotype]])) %>%
pp(p) %>%
param(POSMODE = 2, TSCALE = tscale) %>%
mrgsim(end = years*8766, delta = delta, hmax = 12,
atol = 1e-7, rtol = 1e-6) %>%
as.data.frame()
}
## Freeze the full state vector at a chosen disease year -> the initial
## condition for every fast-timescale experiment below.
msa_state_at <- function(year, phenotype = "MSA_P", p = list(), tscale = 40) {
out <- mod %>%
param(as.list(pheno[[phenotype]])) %>%
pp(p) %>%
param(POSMODE = 2, TSCALE = tscale) %>%
mrgsim(end = year*8766, delta = 24, hmax = 12, atol = 1e-7, rtol = 1e-6) %>%
as.data.frame()
last <- out[nrow(out), ]
cmts <- as.character(outvars(mod)$cmt)
st <- as.list(last[, cmts])
names(st) <- cmts
st
}
## =====================================================================
## STAGE 2a — head-up tilt / active stand from a frozen disease state
## =====================================================================
msa_tilt <- function(state, phenotype = "MSA_P", p = list(),
dose = NULL, settle = 1.0, dur = 0.75) {
m <- mod %>%
param(as.list(pheno[[phenotype]])) %>%
pp(p) %>%
param(POSMODE = 3, TILT_T = settle, MEALAMP = 0, TSCALE = 1, INITMAIN = 0) %>%
init(state)
if (!is.null(dose)) m <- m %>% ev(dose)
m %>% mrgsim(end = settle + dur, delta = 1/240, hmax = 0.002,
atol = 1e-9, rtol = 1e-7, recsort = 3) %>% as.data.frame()
}
## Derive the two diagnostic numbers that define neurogenic OH
tilt_metrics <- function(tt, settle = 1.0) {
base <- tt[which.min(abs(tt$time - (settle - 0.01))), ]
at3 <- tt[which.min(abs(tt$time - (settle + 3/60))), ]
dSBP <- at3$SBP - base$SBP
dHR <- at3$HR - base$HR
data.frame(SBP_supine = base$SBP, SBP_3min = at3$SBP, dSBP = dSBP,
HR_supine = base$HR, HR_3min = at3$HR, dHR = dHR,
HR_SBP_ratio = ifelse(dSBP < 0, dHR/(-dSBP), NA_real_),
nOH = dSBP <= -20,
nOH_strict = dSBP <= -30,
NE_supine = base$NEPL, NE_3min = at3$NEPL,
NE_rise_pct = 100*(at3$NEPL - base$NEPL)/base$NEPL)
}
## =====================================================================
## STAGE 2b — 48-h circadian run with real dosing, from a frozen state
## =====================================================================
msa_day <- function(state, phenotype = "MSA_P", p = list(), dose = NULL,
days = 2, delta = 0.1) {
m <- mod %>%
param(as.list(pheno[[phenotype]])) %>%
pp(p) %>%
param(POSMODE = 2, TSCALE = 1, INITMAIN = 0) %>%
init(state)