Repository navigation
Expand file tree
/
Copy pathptl_mrgsolve_model.R
More file actions
1503 lines (1373 loc) · 78.5 KB
/
Copy pathptl_mrgsolve_model.R
File metadata and controls
1503 lines (1373 loc) · 78.5 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
## Preterm Labor / Spontaneous Preterm Birth (sPTB) — QSP / PK-PD Model
## ============================================================================
## Disease : Spontaneous preterm birth — delivery before 37+0 weeks arising
## not from one disease but from a SYNDROME with several
## aetiologies (intra-amniotic infection/inflammation, functional
## progesterone withdrawal, decidual haemorrhage/thrombin, uterine
## overdistension, cervical insufficiency, an advanced placental
## CRH clock) that all converge on ONE common terminal pathway
## with THREE effector limbs:
## (1) MYOMETRIUM quiescence -> activation -> contraction
## (2) CERVIX softening -> ripening -> dilation
## (3) MEMBRANES weakening -> rupture (PPROM)
##
## THE STRUCTURAL POINT THIS MODEL IS BUILT TO MAKE
## ---------------------------------------------------------------
## Every tocolytic in clinical use acts on limb (1) only. Atosiban
## blocks OXTR, nifedipine blocks Cav1.2, indomethacin blocks COX,
## magnesium antagonises Ca2+ at the myofilament — four different
## molecular targets, all of them inside the same limb, all of them
## leaving limbs (2) and (3) running at full speed. That is not a
## pharmacological accident; it is the reason forty years of
## tocolytic trials have reproducibly bought ~48 hours and
## reproducibly failed to move gestational age at delivery or any
## neonatal outcome. In this model the delivery hazard is written
## as a PRODUCT of the three limbs,
## lambda(t) = LAMBDA0 * f_contr * f_cerv * f_memb * f_term
## so silencing one factor can only reduce lambda by the amount
## that factor contributes. When the cervix is already 12 mm and
## MMP-driven remodelling is running, f_cerv alone keeps lambda
## high and the tocolytic's effect is ARITHMETICALLY BOUNDED —
## PTL_limb_decomposition() prints that bound as a number.
##
## The corollary — and the second thing the model is built to show
## — is that the interventions that DO change outcome are the ones
## that do not treat the uterus at all. Antenatal corticosteroids
## mature the fetal lung; MgSO4 protects the fetal brain; in-utero
## transfer moves the baby to a level-III unit. Tocolysis is
## valuable exactly insofar as it buys the 48 hours those three
## need, and worthless outside that window (hence: no maintenance
## tocolysis). And when the trigger is intra-amniotic infection,
## tocolysis is worse than worthless, because the days are bought
## inside an inflamed amniotic cavity — trading RDS risk for FIRS,
## sepsis and cerebral-palsy risk. PTL_infection_paradox() prices
## this out: prolonging by ~2 days raises cumulative fetal
## inflammatory exposure by ~25%, raises cerebral palsy from 0.060
## to 0.063 and raises early-onset sepsis, while lowering RDS — so
## the composite comes out FLAT. A wash is already a damning
## result for an intervention whose entire justification is that
## the days bought improve outcome.
##
## Authors : Generated by Claude Code Routine (QSP Disease Model Library)
## Engine : mrgsolve (>= 1.5)
##
## Scope : 45 ODE compartments, 186 parameters —
## * Progestogen prevention PK (4): vaginal micronized P4 depot +
## plasma, with an explicit FIRST-UTERINE-PASS local-tissue term
## (plasma concentration is the WRONG exposure metric for this
## route); 17-OHPC IM depot + plasma (flip-flop kinetics).
## * Acute tocolytic PK (8): atosiban 2-compartment IV, nifedipine
## gut + plasma, indomethacin gut + plasma, magnesium central +
## peripheral (renally cleared, CrCl-scaled).
## * Fetal-directed drug PK (7): betamethasone phosphate (fast
## ester) and acetate (slow depot) IM, maternal plasma, fetal
## plasma across an 11b-HSD2 barrier betamethasone evades;
## terbutaline depot + plasma + a beta2-desensitisation state.
## * Endocrine ignition (4): placental CRH, fetal DHEA-S, fetal
## cortisol, estriol.
## * Progesterone-withdrawal axis (1): PR-A/PR-B index.
## * Infection / inflammation (4): intra-amniotic bacterial load,
## NF-kB activity, amniotic IL-6, fetal inflammatory response.
## * Uterotonin machinery (5): COX-2, prostaglandin pool, OXTR
## density, connexin-43, myometrial cAMP.
## * Contraction (2): intracellular Ca2+ index, contraction index.
## * Cervix (4): cross-linked collagen, hyaluronan/hydration,
## MMP activity, sonographic cervical length.
## * Membranes (2): tensile strength, cervicovaginal fFN.
## * Fetal maturation (2): surfactant pool, brain Mg exposure.
## * Outcome accounting (2): cumulative delivery hazard,
## cumulative fetal inflammatory exposure.
##
## Time : DAYS. t = 0 is a user-set gestational age (GA0_WK; 20 weeks for
## prevention scenarios, 27-30 weeks for acute scenarios).
## Gestational age in weeks is GAW = GA0_WK + t/7.
## Default horizons: 140 d for prevention (20 -> 40 wk),
## 14-21 d for acute threatened preterm labour.
##
## Units : ng/mL (progesterone, 17-OHPC, atosiban, nifedipine,
## indomethacin, betamethasone, terbutaline, fFN),
## mg/dL (serum magnesium), pg/mL (CRH), ng/mL (amniotic IL-6),
## mm (cervical length), log10 CFU/mL (intra-amniotic bacteria),
## MVU-like index (contraction), and fold-of-baseline or 0-1
## normalised indices for receptor densities, transcription
## factors and tissue-integrity states. Those normalised states
## are relative-scale QSP variables, not assayable concentrations
## — this library's standard convention where the biology is well
## described qualitatively but not calibrated absolutely in vivo.
##
## Scenarios: 12 prebuilt — see PTL_scenarios() / PTL_simulate_scenarios().
## Analyses : PTL_ledger(), PTL_limb_decomposition(), PTL_acs_window(),
## PTL_infection_paradox(), PTL_tocolytic_head_to_head().
##
## Anchors :
## * McLean 1995 Nat Med 1:460 — placental CRH rises exponentially with a
## doubling time ~3.3 wk and the trajectory is set weeks to months before
## delivery ("the placental clock"). Sets CRH_TDOUBLE and CLOCK_ADV.
## * Challis 2000 Endocr Rev 21:514; Smith 2007 Sci Am / 2009 — CRH ->
## fetal adrenal DHEA-S -> 16a-OH-DHEAS -> placental aromatase ->
## estriol; fetal cortisol feeds back POSITIVELY on placental CRH
## (opposite in sign to the hypothalamic loop). Sets FB_CORT.
## * Mesiano 2002 JCEM 87:2924; Merlino 2007 JCEM 92:1927 — functional
## progesterone withdrawal is a rise in the myometrial PR-A/PR-B ratio,
## NOT a fall in plasma progesterone. This is why PRW is a state
## variable rather than a function of plasma P4.
## * Renthal 2010 PNAS 107:20828 — miR-200 / ZEB1-2 switch releasing the
## CAP genes; motivates the sigmoidal PRW -> CAP coupling.
## * Romero 2014 Am J Reprod Immunol 72:458; Romero 2015 J Perinat Med
## 43:19 — MIAC vs STERILE intra-amniotic inflammation (the latter about
## twice as common); amniotic-fluid IL-6 > 2.6 ng/mL threshold.
## * Gomez 1998 Am J Obstet Gynecol 179:194 — FIRS (cord IL-6 > 11 pg/mL)
## and its independent association with neonatal morbidity.
## * Iams 1996 NEJM 334:567 — mid-trimester cervical-length distribution
## (median 35 mm, 10th centile 25 mm) and the continuous, inverse,
## non-linear CL -> PTB risk relation. Sets CL0 and the f_cerv form.
## * Hassan 2011 Ultrasound Obstet Gynecol 38:18 (PREGNANT) — vaginal P4
## in CL 10-20 mm: birth < 33 wk 8.9% vs 16.1%, RR 0.55.
## * EPPPIC Collaborative 2021 Lancet 397:1183 — IPD meta-analysis:
## vaginal progesterone reduces birth < 34 wk in singletons at risk
## (RR ~0.78); no clear benefit in unselected twins. Calibration target
## for EMAX_P4_PRW / EC50_P4_LOC / FUP_RATIO.
## * Meis 2003 NEJM 348:2379 (positive) vs Blackwell 2020 Am J Perinatol
## 37:127 (PROLONG, null) — 17-OHPC. Encoded through the SAME PD term as
## vaginal P4 but with a much lower effective local exposure
## (OHPC_LOCAL), which is one mechanistic reading of the discrepancy.
## The model does NOT claim to resolve what the trials did not.
## * Goodwin 1995 Am J Obstet Gynecol 173:913; Worldwide Atosiban vs
## Beta-agonists Study Group 2001 BJOG 108:133 — atosiban PK: Vc ~18 L,
## CL ~42 L/h, t1/2alpha ~13 min, t1/2beta ~1.7 h; 6.75 mg bolus,
## 300 ug/min x 3 h, then 100 ug/min for up to 45 h.
## * Flenady 2014 Cochrane CD002255 (calcium-channel blockers) and
## CD004452 (oxytocin-receptor antagonists) — tocolytics delay birth by
## 48 h / 7 d without improving neonatal outcome. The multiplicative
## hazard exists to reproduce exactly this dissociation.
## * Roberts 2020 Cochrane CD004454 — antenatal corticosteroids: RDS
## RR 0.66, neonatal death RR 0.78, severe IVH RR 0.55.
## * Gyamfi-Bannerman 2016 NEJM 374:1311 (ALPS) — late-preterm steroids:
## respiratory morbidity 11.6% vs 14.4%; neonatal hypoglycaemia 24% vs
## 15%. The hypoglycaemia trade-off is exported as a readout.
## * Rouse 2008 NEJM 359:895 (BEAM); Doyle 2009 Cochrane CD004661 — MgSO4
## neuroprotection < 32 wk: moderate/severe CP 1.9% vs 3.5%, RR ~0.68,
## NNT ~63. Sets EMAX_MGN / EC50_MGN.
## * Kenyon 2001 Lancet 357:979 (ORACLE-I) — erythromycin after PPROM
## prolongs latency; co-amoxiclav increases neonatal NEC.
## * Moise 1988 NEJM 319:327 — indomethacin and fetal ductal constriction,
## steeply gestational-age dependent. Sets DUCT_GA50 = 31.5 wk.
## * Stoll 2015 JAMA 314:1039; Manuck 2016 Am J Obstet Gynecol 215:103.e1 —
## gestational-age-specific neonatal survival and morbidity, used for the
## RDS / IVH / NEC / death logistics.
##
## Disclaimer: Research / education / hypothesis generation only. This is NOT
## a tocolysis protocol, NOT a dosing calculator, and NOT validated
## against patient data. Threatened preterm labour is an obstetric
## emergency requiring clinical judgment.
## ============================================================================
[PROB]
QSP model of spontaneous preterm birth. Four aetiologic ignition arms (placental
CRH clock, functional progesterone withdrawal, intra-amniotic infection/sterile
inflammation, decidual thrombin and uterine overdistension) converge on a common
terminal pathway whose three effector limbs — myometrial contraction, cervical
ripening and membrane weakening — enter the delivery hazard MULTIPLICATIVELY.
PK/PD is included for vaginal progesterone, 17-OHPC, atosiban, nifedipine,
indomethacin, magnesium sulfate, terbutaline, betamethasone and erythromycin,
with gestational-age-dependent neonatal outcome models for RDS, severe IVH,
NEC, early-onset sepsis, cerebral palsy and death.
[PLUGIN] Rcpp
[GLOBAL]
// Variables computed once per record in $MAIN and reused in $ODE / $TABLE.
// (Locals declared inside $MAIN are not visible to $ODE, so they live here.)
double CLMG; // magnesium clearance, CrCl-scaled (L/day)
double STRETCH; // total uterine overdistension load (0-1+)
double PRWINIT; // starting functional-withdrawal index
[PARAM] @annotated
//--- Patient / pregnancy profile -------------------------------------------
GA0_WK : 20.0 : weeks, gestational age at t = 0
WT : 75.0 : kg, maternal weight in pregnancy
CRCL : 130.0 : mL/min, maternal creatinine clearance (gestational hyperfiltration)
CL0 : 35.0 : mm, cervical length at t = 0
NFETUS : 1.0 : number of fetuses (1 = singleton, 2 = twins)
PRIOR_PTB : 0.0 : 0/1, prior spontaneous preterm birth
PPROM_FLAG : 0.0 : 0/1, membranes already ruptured at t = 0
ACTIVATED : 0.0 : 0-1, myometrial activation ALREADY present at t = 0
(a woman presenting in preterm labour does not start
from quiescence — set ~0.7 for acute scenarios)
//--- Aetiologic driver strengths (the four ignition arms) -------------------
CLOCK_ADV : 0.0 : weeks, advancement of the placental CRH clock set-point
INFECT_DRIVE : 0.0 : 0-1, intra-amniotic microbial inoculum intensity
STERILE_DRIVE: 0.05: 0-1, sterile (DAMP-driven) inflammatory drive
THROMBIN_DRV : 0.0 : 0-1, decidual haemorrhage / thrombin arm intensity
STRETCH_DRV : 0.0 : 0-1, overdistension beyond that implied by NFETUS
CERV_FRAGILE : 0.0 : 0-1, intrinsic cervical insufficiency (LEEP, anomaly)
//--- Placental CRH clock ----------------------------------------------------
CRHBASE : 60.0 : pg/mL, plasma CRH at 20 weeks
CRH_TDOUBLE : 23.1 : days, CRH doubling time (3.3 weeks; McLean 1995)
CRH_MAX : 3000.0: pg/mL, CRH ceiling near term
KOUT_CRH : 0.35: 1/day, CRH elimination scalar
FB_CORT : 0.55: -, strength of fetal-cortisol -> placental-CRH POSITIVE feedback
//--- Fetal adrenal / estriol axis -------------------------------------------
KIN_DHEAS : 0.9 : 1/day, CRH-driven fetal DHEA-S production scalar
KOUT_DHEAS : 0.9 : 1/day, DHEA-S turnover
EC50_CRH_ADR : 600.0: pg/mL, CRH potency at the fetal adrenal
KIN_CORTF : 0.8 : 1/day, fetal cortisol production scalar
KOUT_CORTF : 0.8 : 1/day, fetal cortisol turnover
KIN_E3 : 1.2 : 1/day, DHEA-S -> estriol aromatisation scalar
KOUT_E3 : 1.2 : 1/day, estriol turnover
//--- Functional progesterone withdrawal (PR-A / PR-B) -----------------------
PRWBASE : 0.15: -, baseline functional-withdrawal index (0-1)
KIN_PRW : 0.055: 1/day, drive toward withdrawal
KOUT_PRW : 0.055: 1/day, reversal of withdrawal
E3_ON_PRW : 0.45: -, estriol contribution to withdrawal
INFL_ON_PRW : 0.60: -, inflammatory contribution (IL-1/TNF -> 20a-HSD, PR-A)
EMAX_P4_PRW : 0.55: -, maximal suppression of withdrawal by local progestogen
EC50_P4_LOC : 120.0 : ng/mL-eq, local (endometrial) progestogen potency
P4_NFKB : 0.60: -, maximal PR-mediated suppression of the NF-kB drive
P4_CERV : 0.55: -, maximal progestogen suppression of cervical MMP induction
//--- Intra-amniotic infection & inflammation --------------------------------
BACTBASE : 0.0 : log10 CFU/mL, intra-amniotic bacterial load at t = 0
KGROW_BACT : 0.55: 1/day, net bacterial growth rate in amniotic fluid
BACT_MAX : 7.0 : log10 CFU/mL, carrying capacity
KKILL_ABX : 0.60: 1/day per unit antibiotic effect (partial — Ureaplasma tolerance)
IL6BASE : 0.4 : ng/mL, baseline amniotic-fluid IL-6
KIN_IL6 : 12.0 : ng/mL/day, maximal NF-kB-driven IL-6 production
KOUT_IL6 : 1.4 : 1/day, IL-6 clearance from amniotic fluid
IL6_THRESH : 2.6 : ng/mL, intra-amniotic inflammation threshold (reporting only)
KIN_NFKB : 3.0 : 1/day, NF-kB activation
KOUT_NFKB : 3.0 : 1/day, NF-kB deactivation
EC50_BACT : 3.5 : log10 CFU/mL, PAMP potency at TLR2/4
PR_ON_NFKB : 0.50: -, PR-mediated trans-repression of NF-kB (lost with withdrawal)
KIN_FIRS : 0.18: 1/day, fetal inflammatory response induction
KOUT_FIRS : 0.05: 1/day, FIRS resolution — deliberately SLOW, so that
fetal inflammatory injury ACCUMULATES with exposure
time. This is what makes prolonging a pregnancy
inside an infected uterus cost something.
EC50_IL6_FIRS: 4.0 : ng/mL, amniotic IL-6 potency for FIRS
//--- Uterotonin machinery ---------------------------------------------------
COX2BASE : 1.0 : fold, basal myometrial/amnion COX-2
KIN_COX2 : 1.8 : 1/day, COX-2 induction scalar
KOUT_COX2 : 1.8 : 1/day, COX-2 turnover
COX2_MAX : 12.0 : fold, maximal COX-2 induction
PGBASE : 1.0 : fold, basal prostaglandin pool
KIN_PG : 2.4 : 1/day, PG synthesis scalar
KOUT_PG : 2.4 : 1/day, PG disposal scalar
PGDH_FRAC : 0.70: -, fraction of PG disposal via the chorionic 15-PGDH barrier
KM_PG : 8.0 : fold, PG pool giving half-maximal uterotonic signalling
OXTRBASE : 1.0 : fold, basal OXTR density
KIN_OXTR : 0.30: 1/day, OXTR upregulation scalar
KOUT_OXTR : 0.30: 1/day, OXTR turnover
OXTR_MAX : 120.0 : fold, maximal term OXTR upregulation (reported 100-200x)
KM_OXTR : 25.0 : fold, OXTR density giving half-maximal oxytocin signalling
CX43BASE : 1.0 : fold, basal connexin-43
KIN_CX43 : 0.45: 1/day, Cx43 induction scalar
KOUT_CX43 : 0.45: 1/day, Cx43 turnover
CX43_MAX : 18.0 : fold, maximal Cx43 induction
OT_TONE : 1.0 : fold, endogenous oxytocin tone multiplier
//--- Myometrial excitation-contraction coupling ------------------------------
CAIBASE : 1.0 : fold, basal myometrial Ca2+ index
KOUT_CAI : 60.0 : 1/day, Ca2+ index equilibration (fast on the day scale)
CAMPBASE : 1.0 : fold, basal myometrial cAMP
KOUT_CAMP : 90.0 : 1/day, cAMP equilibration (fast)
KIN_CONTR : 36.0 : 1/day, contraction-index equilibration
CONTR_MAX : 260.0 : MVU-like units, maximal contraction index
HILL_CONTR : 2.6 : -, Hill coefficient of the Ca2+ -> contraction relation
EC50_CONTR : 3.0 : fold, Ca2+ index giving half-maximal contraction
W_TONIC : 0.30: -, tonic (agonist-independent) component of the Ca2+ drive
W_OT : 1.80: -, weight of the oxytocin arm in the Ca2+ drive
W_PG : 3.00: -, weight of the prostaglandin arm in the Ca2+ drive
W_PAR1 : 1.20: -, weight of the thrombin/PAR-1 arm in the Ca2+ drive
CAMP_BRAKE : 0.80: -, maximal fractional suppression of contraction by cAMP
EC50_CAMP : 2.5 : fold, cAMP potency as a contraction brake
//--- Cervical remodelling ----------------------------------------------------
COLLBASE : 1.0 : fold, cross-linked collagen integrity at t = 0
KDEG_COLL : 0.09: 1/day, maximal MMP-driven collagen degradation
HILL_MMPCOL : 2.0 : -, cooperativity of MMP-driven collagen loss (a threshold effect)
KREP_COLL : 0.020: 1/day, collagen repair / synthesis
HABASE : 0.10: -, cervical hyaluronan / hydration index at t = 0
KIN_HA : 0.30: 1/day, HAS2-driven hyaluronan accumulation
KOUT_HA : 0.30: 1/day, hyaluronan turnover
MMPBASE : 1.0 : fold, cervical & membrane MMP activity at t = 0
KIN_MMP : 1.5 : 1/day, MMP induction
KOUT_MMP : 1.5 : 1/day, MMP turnover / TIMP restraint
MMP_MAX : 14.0 : fold, maximal MMP induction
CL_MIN : 3.0 : mm, functional floor of cervical length
KCL : 0.22: 1/day, rate at which CL tracks its structural target
CERC_SUPPORT : 0.0 : 0/1, cervical cerclage in place
CERC_EFF : 7.0 : mm, effective CL added by a well-placed cerclage
PESS_SUPPORT : 0.0 : 0/1, Arabin pessary in place
PESS_EFF : 3.0 : mm, effective CL added by a pessary
//--- Fetal membranes & fFN ---------------------------------------------------
MEMBBASE : 1.0 : fold, membrane tensile strength in intact membranes
KDEG_MEMB : 0.11: 1/day, maximal MMP-9 / apoptosis-driven weakening
KREP_MEMB : 0.010: 1/day, membrane repair (minimal — membranes barely heal)
FFNBASE : 8.0 : ng/mL, baseline cervicovaginal fetal fibronectin
KIN_FFN : 260.0 : ng/mL/day, fFN release scalar with choriodecidual disruption
KOUT_FFN : 1.1 : 1/day, fFN clearance from the cervicovaginal pool
FFN_POS : 50.0 : ng/mL, standard positivity threshold (reporting only)
//--- Delivery hazard ----------------------------------------------------------
LAMBDA0 : 0.0024: 1/day, baseline hazard scalar
CONTR_FLOOR : 0.15: -, irreducible floor of the contraction hazard limb
CONTR_REF : 40.0 : MVU-like units, contraction index giving a unit hazard factor
CERV_FLOOR : 0.15: -, floor of the cervical hazard limb (a long closed
cervix is genuinely PROTECTIVE, not merely neutral)
CERV_AMP : 14.0 : -, amplitude of the cervical hazard limb
CL_REF : 14.0 : mm, cervical length at half-maximal cervical hazard
CL_STEEP : 3.0 : mm, steepness of the cervical hazard term
MEMB_AMP : 4.0 : -, amplitude of the membrane hazard limb
MEMB_REF : 0.45: fold, membrane strength at half-maximal membrane hazard
MEMB_STEEP : 0.10: fold, steepness of the membrane hazard term
GA_RIPEN : 38.0 : weeks, GA at which the physiologic term drive is half-maximal
GA_STEEP : 0.8 : weeks, steepness of the physiologic term drive
TERM_AMP : 150.0 : -, amplitude of the physiologic term hazard multiplier
PPROM_HAZ : 9.0 : -, hazard multiplier once membranes have ruptured
//--- Fetal lung maturation ----------------------------------------------------
SURFBASE : 0.05: -, surfactant / lung-maturity index at 20 weeks
KIN_SURF : 0.055: 1/day, endogenous (cortisol-driven) maturation rate
GA_SURF50 : 34.5 : weeks, GA at half-maximal endogenous surfactant
GA_SURF_ST : 2.2 : weeks, steepness of the endogenous maturation curve
EMAX_SURF : 0.55: -, maximal betamethasone-driven increment in lung maturity
EC50_BETF : 6.0 : ng/mL, fetal betamethasone potency at the lung GR
KSURF_ACS : 0.9 : 1/day, rate of the corticosteroid-driven surfactant response
ACS_DECAY : 0.10: 1/day, waning of the corticosteroid effect
//--- Drug: vaginal micronized progesterone ------------------------------------
KA_P4V : 3.5 : 1/day, vaginal absorption rate
F_P4V : 0.03: -, systemic bioavailability of vaginal P4 (low by design)
VD_P4 : 120.0 : L, apparent volume of distribution
CL_P4 : 480.0 : L/day, progesterone clearance
FUP_RATIO : 28.0 : -, first-uterine-pass local:plasma concentration ratio
//--- Drug: 17-alpha-hydroxyprogesterone caproate -------------------------------
KA_OHPC : 0.05: 1/day, IM depot absorption (flip-flop; absorption-limited)
F_OHPC : 0.85: -, IM bioavailability
VD_OHPC : 500.0 : L, apparent volume of distribution
CL_OHPC : 2000.0: L/day, apparent clearance
OHPC_LOCAL : 6.0 : -, myometrial:plasma partitioning (no first-uterine pass)
//--- Drug: atosiban (2-compartment IV) ------------------------------------------
VC_ATO : 18.0 : L, central volume
VP_ATO : 12.0 : L, peripheral volume
CL_ATO : 1010.0: L/day, systemic clearance (~42 L/h)
Q_ATO : 360.0 : L/day, intercompartmental clearance
KI_ATO : 90.0 : ng/mL, competitive OXTR inhibition constant
//--- Drug: nifedipine ------------------------------------------------------------
KA_NIF : 14.0 : 1/day, oral absorption
F_NIF : 0.50: -, oral bioavailability (extensive first pass)
VD_NIF : 55.0 : L, apparent volume of distribution
CL_NIF : 620.0 : L/day, clearance (CYP3A4)
IC50_NIF : 35.0 : ng/mL, Cav1.2 blockade potency
IMAX_NIF : 0.72: -, maximal fractional block of Ca2+ influx
//--- Drug: indomethacin -----------------------------------------------------------
KA_IND : 16.0 : 1/day, rectal / oral absorption
F_IND : 0.90: -, bioavailability
VD_IND : 38.0 : L, apparent volume of distribution
CL_IND : 115.0 : L/day, clearance
IC50_IND : 450.0 : ng/mL, COX inhibition potency
IMAX_IND : 0.85: -, maximal fractional COX inhibition
F_PLAC_IND : 0.50: -, fetal:maternal indomethacin ratio
DUCT_GA50 : 31.5 : weeks, GA at half-maximal ductal-constriction susceptibility
DUCT_GA_ST : 1.2 : weeks, steepness of the ductal susceptibility curve
//--- Drug: magnesium sulfate --------------------------------------------------------
MG_BASE : 1.9 : mg/dL, baseline serum magnesium
VC_MG : 16.0 : L, central (extracellular) volume
VP_MG : 20.0 : L, peripheral volume
Q_MG : 30.0 : L/day, intercompartmental clearance
CL_MG_REF : 95.0 : L/day, magnesium clearance at CrCl 130 mL/min (renal only)
IC50_MG : 5.5 : mg/dL, Mg potency as a Ca2+ / MLCK antagonist
IMAX_MG : 0.45: -, maximal fractional contraction suppression by Mg
EC50_MGN : 1.2 : mg/dL ABOVE BASELINE, Mg potency for fetal neuroprotection
KMG_BRAIN : 1.8 : 1/day, equilibration of fetal brain Mg exposure
//--- Drug: betamethasone -------------------------------------------------------------
KA_BETP : 22.0 : 1/day, betamethasone SODIUM PHOSPHATE absorption (fast ester)
KA_BETA : 1.6 : 1/day, betamethasone ACETATE absorption (slow depot)
VD_BET : 60.0 : L, apparent volume of distribution
CL_BET : 185.0 : L/day, maternal clearance (t1/2 ~5.5 h)
F_PLAC_BET : 0.30: -, fetal:maternal ratio (betamethasone escapes 11b-HSD2)
KEQ_BETF : 12.0 : 1/day, fetal compartment equilibration
//--- Drug: terbutaline (beta2 agonist; largely historical) -----------------------------
KA_TERB : 30.0 : 1/day, subcutaneous absorption
VD_TERB : 110.0 : L, apparent volume of distribution
CL_TERB : 200.0 : L/day, clearance
EC50_TERB : 3.5 : ng/mL, beta2 potency
EMAX_TERB : 3.0 : fold, maximal cAMP elevation
KDES_ON : 1.1 : 1/day, beta2 desensitisation development
KDES_OFF : 0.30: 1/day, resensitisation
DES_MAX : 0.85: -, maximal loss of beta2 responsiveness (tachyphylaxis)
//--- Adjuncts (switches) --------------------------------------------------------------
ERY_ON : 0.0 : 0/1, erythromycin 250 mg qid for PPROM (ORACLE-I)
COAMOX_ON : 0.0 : 0/1, co-amoxiclav (records the NEC penalty; NOT recommended)
AZI_ON : 0.0 : 0/1, azithromycin anti-Ureaplasma (investigational)
TRANSFER_ON : 0.0 : 0/1, in-utero transfer to a level-III centre
[CMT] @annotated
//--- Progestogen PK --------------------------------------------------------
P4V_DEP : vaginal micronized progesterone depot (mg)
P4_C : plasma progesterone amount, exogenous only (ng/mL x L)
OHPC_DEP : 17-OHPC IM depot (mg)
OHPC_C : 17-OHPC plasma amount (mg)
//--- Acute tocolytic PK ----------------------------------------------------
ATO_C : atosiban central amount (mg)
ATO_P : atosiban peripheral amount (mg)
NIF_GUT : nifedipine gut amount (mg)
NIF_C : nifedipine plasma amount (mg)
IND_GUT : indomethacin gut / rectal amount (mg)
IND_C : indomethacin plasma amount (mg)
MG_C : magnesium central amount above baseline (mg/dL x L)
MG_P : magnesium peripheral amount (mg/dL x L)
//--- Fetal-directed drug PK ------------------------------------------------
BET_DP : betamethasone phosphate ester depot (mg)
BET_DA : betamethasone acetate slow depot (mg)
BET_C : betamethasone maternal plasma amount (mg)
BET_F : betamethasone fetal plasma concentration (ng/mL)
TERB_DEP : terbutaline subcutaneous depot (mg)
TERB_C : terbutaline plasma amount (mg)
DESENS : beta2 receptor desensitisation index (0-1)
//--- Endocrine ignition ----------------------------------------------------
CRH : placental CRH (pg/mL)
DHEAS : fetal adrenal DHEA-S (fold of 20-week baseline)
CORTF : fetal cortisol (fold of 20-week baseline)
E3 : estriol (fold of 20-week baseline)
//--- Progesterone withdrawal -----------------------------------------------
PRW : PR-A/PR-B functional progesterone withdrawal index (0-1)
//--- Infection / inflammation ----------------------------------------------
BACT : intra-amniotic bacterial load (log10 CFU/mL)
NFKB : myometrial / decidual NF-kB activity (0-1)
IL6 : amniotic fluid IL-6 (ng/mL)
FIRS : fetal inflammatory response index (0-1)
//--- Uterotonin machinery --------------------------------------------------
COX2 : COX-2 expression (fold)
PG : prostaglandin pool (fold)
OXTR : oxytocin receptor density (fold)
CX43 : connexin-43 gap junction density (fold)
CAMP : myometrial cAMP (fold)
//--- Contraction -----------------------------------------------------------
CAI : myometrial intracellular Ca2+ index (fold)
CONTR : contraction index (MVU-like)
//--- Cervix ----------------------------------------------------------------
COLL : cervical cross-linked collagen integrity (fold)
HA : cervical hyaluronan / hydration index (0-1)
MMP : cervical & membrane MMP activity (fold)
CLEN : sonographic cervical length (mm)
//--- Membranes -------------------------------------------------------------
MEMB : fetal membrane tensile strength (fold)
FFN : cervicovaginal fetal fibronectin (ng/mL)
//--- Fetal maturation ------------------------------------------------------
SURF : fetal lung surfactant / maturity index (0-1)
MGBRAIN : fetal brain magnesium exposure index (mg/dL-equivalent)
//--- Outcome accounting ----------------------------------------------------
CUMHAZ : cumulative delivery hazard (dimensionless)
CUMINFL : cumulative fetal inflammatory exposure (ng/mL x day)
[MAIN]
// ---- individual scaling ---------------------------------------------------
CLMG = CL_MG_REF * (CRCL / 130.0); // magnesium is cleared renally, and only renally
STRETCH = STRETCH_DRV + 0.55 * (NFETUS - 1.0); // twins carry an intrinsic stretch load
PRWINIT = PRWBASE + 0.05 * PRIOR_PTB + 0.10 * CERV_FRAGILE;
// ---- initial conditions ---------------------------------------------------
CRH_0 = CRHBASE;
DHEAS_0 = 1.0;
CORTF_0 = 1.0;
E3_0 = 1.0;
// A woman presenting in threatened preterm labour is not starting from
// quiescence: ACTIVATED blends the initial state toward an activated
// myometrium. If the underlying drive is weak the state relaxes back — which
// is exactly why most women with threatened preterm labour do not deliver.
PRW_0 = PRWINIT + ACTIVATED * (0.85 - PRWINIT);
BACT_0 = BACTBASE;
NFKB_0 = 0.05 + 0.50 * STERILE_DRIVE + 0.25 * ACTIVATED;
IL6_0 = IL6BASE;
FIRS_0 = 0.0;
COX2_0 = COX2BASE;
PG_0 = PGBASE;
OXTR_0 = OXTRBASE + ACTIVATED * (0.45 * OXTR_MAX - OXTRBASE);
CX43_0 = CX43BASE + ACTIVATED * (0.55 * CX43_MAX - CX43BASE);
CAMP_0 = CAMPBASE;
CAI_0 = CAIBASE;
CONTR_0 = 4.0 + 80.0 * ACTIVATED;
COLL_0 = COLLBASE;
HA_0 = HABASE;
MMP_0 = MMPBASE + ACTIVATED * 0.25 * MMP_MAX;
CLEN_0 = CL0;
MEMB_0 = MEMBBASE * (1.0 - 0.55 * PPROM_FLAG);
FFN_0 = FFNBASE;
SURF_0 = SURFBASE;
DESENS_0 = 0.0;
MGBRAIN_0 = 0.0;
CUMHAZ_0 = 0.0;
CUMINFL_0 = 0.0;
[ODE]
// =========================================================================
// 0. Gestational age and the physiologic term drive
// =========================================================================
double GA = GA0_WK + SOLVERTIME / 7.0;
double TERMDRIVE = 1.0 / (1.0 + exp(-(GA - GA_RIPEN) / GA_STEEP));
// =========================================================================
// 1. Drug concentrations derived from amounts
// =========================================================================
double CP4 = P4_C / VD_P4; // ng/mL, exogenous progesterone
double CP4_LOC = CP4 * FUP_RATIO; // first-uterine-pass local exposure
double COHPC = 1000.0 * OHPC_C / VD_OHPC; // ng/mL 17-OHPC
double CP4_EFF = CP4_LOC + OHPC_LOCAL * COHPC; // combined LOCAL progestogen exposure
double CATO = 1000.0 * ATO_C / VC_ATO; // ng/mL atosiban
double CNIF = 1000.0 * NIF_C / VD_NIF; // ng/mL nifedipine
double CIND = 1000.0 * IND_C / VD_IND; // ng/mL indomethacin
double CMG = MG_BASE + MG_C / VC_MG; // mg/dL serum magnesium
double CBET = 1000.0 * BET_C / VD_BET; // ng/mL maternal betamethasone
double CTERB = 1000.0 * TERB_C / VD_TERB; // ng/mL terbutaline
// Local progestogen protection, computed once — it acts on THREE things:
// the PR-A/PR-B withdrawal drive, the PR/NF-kB trans-repression, and cervical
// MMP induction. That third arm is why progestogens are a CERVICAL therapy as
// much as a myometrial one, and why they work in short cervix specifically.
double p4prot = EMAX_P4_PRW * CP4_EFF / (EC50_P4_LOC + CP4_EFF);
// =========================================================================
// 2. ARM 1 — the placental CRH clock
// Exponential growth to a ceiling, accelerated by a user-set clock
// advancement, by ischaemia/stress proxies, and by the fetal cortisol
// POSITIVE feedback that is unique to the placenta.
// =========================================================================
double kgrow = log(2.0) / CRH_TDOUBLE;
double advance = exp(0.16 * CLOCK_ADV);
double crhfb = 1.0 + FB_CORT * (CORTF - 1.0);
dxdt_CRH = kgrow * advance * crhfb * CRH * (1.0 - CRH / CRH_MAX)
+ 0.10 * KOUT_CRH * CRHBASE * (1.0 + 2.0 * THROMBIN_DRV)
- 0.10 * KOUT_CRH * CRH;
double crhdrive = CRH / (EC50_CRH_ADR + CRH);
dxdt_DHEAS = KIN_DHEAS * (0.25 + 3.5 * crhdrive) - KOUT_DHEAS * DHEAS;
dxdt_CORTF = KIN_CORTF * (0.30 + 2.6 * crhdrive) - KOUT_CORTF * CORTF;
dxdt_E3 = KIN_E3 * DHEAS - KOUT_E3 * E3;
// =========================================================================
// 3. ARM 3 — intra-amniotic infection and inflammation
// =========================================================================
// Antibiotic effect is deliberately PARTIAL: ORACLE-I bought latency, not
// sterilisation, and Ureaplasma is only intermittently eradicated in vivo.
double abx = KKILL_ABX * (ERY_ON + 0.8 * COAMOX_ON + 1.4 * AZI_ON);
// Ruptured membranes remove the barrier to ascending colonisation.
double inoc = INFECT_DRIVE * (1.0 + 0.8 * (MEMB < 0.35 ? 1.0 : 0.0));
dxdt_BACT = KGROW_BACT * inoc * (1.0 - BACT / BACT_MAX) * (BACT + 0.35 * inoc)
- abx * BACT;
double pamp = BACT / (EC50_BACT + BACT);
double damp = STERILE_DRIVE + 0.45 * THROMBIN_DRV + 0.30 * STRETCH
+ 0.25 * (1.0 - MEMB);
// Intact PR-B trans-represses NF-kB; functional withdrawal removes that brake.
double prrepress = PR_ON_NFKB * (1.0 - PRW);
double nfdrive = (pamp + damp) * (1.0 - prrepress) * (1.0 - P4_NFKB * p4prot);
dxdt_NFKB = KIN_NFKB * nfdrive * (1.0 - NFKB) - KOUT_NFKB * NFKB;
dxdt_IL6 = KIN_IL6 * NFKB - KOUT_IL6 * IL6;
double firsdrive = IL6 / (EC50_IL6_FIRS + IL6);
dxdt_FIRS = KIN_FIRS * firsdrive * (1.0 - FIRS) - KOUT_FIRS * FIRS;
// =========================================================================
// 4. ARM 2 — functional progesterone withdrawal
// Withdrawal is a PR-A/PR-B RATIO shift, not a fall in plasma P4.
// Progestogen therapy works by holding PR-B signalling, which is why the
// LOCAL (first-uterine-pass) exposure is the concentration that matters.
// =========================================================================
double e3term = (E3 - 1.0) / (1.0 + (E3 > 1.0 ? E3 - 1.0 : 0.0));
if (e3term < 0.0) e3term = 0.0;
double prwdrive = E3_ON_PRW * e3term + INFL_ON_PRW * NFKB
+ 0.35 * TERMDRIVE + 0.20 * STRETCH;
dxdt_PRW = KIN_PRW * prwdrive * (1.0 - p4prot) * (1.0 - PRW)
- KOUT_PRW * PRW * (0.4 + p4prot);
// =========================================================================
// 5. The CAP programme — uterotonin machinery
// =========================================================================
double capdrive = 0.55 * PRW + 0.35 * NFKB + 0.30 * TERMDRIVE
+ 0.25 * STRETCH + 0.20 * THROMBIN_DRV;
// COX-2 EXPRESSION is NF-kB / CAP driven ...
dxdt_COX2 = KIN_COX2 * (0.6 * NFKB + 0.5 * capdrive) * (COX2_MAX - COX2)
- KOUT_COX2 * COX2;
// ... while indomethacin inhibits COX ACTIVITY, not expression.
double coxinhib = IMAX_IND * CIND / (IC50_IND + CIND);
double coxact = COX2 * (1.0 - coxinhib);
// Progestogen maintains the chorionic 15-PGDH disposal barrier, so withdrawal
// simultaneously raises PG synthesis AND removes PG disposal — a double hit.
double pgdh = PGDH_FRAC * (1.0 - 0.75 * PRW) * (1.0 - 0.5 * NFKB);
dxdt_PG = KIN_PG * 8.0 * coxact / COX2_MAX
- KOUT_PG * PG * ((1.0 - PGDH_FRAC) + pgdh);
// OXTR and Cx43 — the two CAPs that convert an activated myometrium into a
// synchronised, oxytocin-responsive one.
dxdt_OXTR = KIN_OXTR * capdrive * (OXTR_MAX - OXTR) - KOUT_OXTR * OXTR;
dxdt_CX43 = KIN_CX43 * (capdrive + 0.3 * PG / (1.0 + PG)) * (CX43_MAX - CX43)
- KOUT_CX43 * CX43;
// =========================================================================
// 6. Relaxation: cAMP (the beta2-agonist target, with tachyphylaxis)
// =========================================================================
double b2resp = 1.0 - DESENS;
double terbocc = CTERB / (EC50_TERB + CTERB);
double terbeff = EMAX_TERB * b2resp * terbocc;
dxdt_DESENS = KDES_ON * terbocc * (DES_MAX - DESENS) - KDES_OFF * DESENS;
dxdt_CAMP = KOUT_CAMP * (CAMPBASE + terbeff - CAMP);
// =========================================================================
// 7. Excitation-contraction coupling
// Ca2+ is driven by oxytocin (via OXTR), prostaglandin (FP/EP1), thrombin
// (PAR-1) and electrical coupling (Cx43); it is opposed by nifedipine
// (Cav1.2 block) and magnesium (competitive Ca2+ antagonism).
// Atosiban enters as a COMPETITIVE OXTR antagonist — 1/(1 + C/Ki) — not
// as an Emax term, because that is the form that reproduces the clinical
// observation that atosiban's apparent efficacy falls as endogenous
// oxytocin tone rises. An Emax form cannot express that.
// =========================================================================
double atoblock = 1.0 / (1.0 + CATO / KI_ATO);
// Receptor DENSITY is not signal: OXTR rises 100x but the downstream Gq
// response saturates, so the oxytocin arm enters as a fractional occupancy.
double ot_sig = OT_TONE * (OXTR / (KM_OXTR + OXTR)) * atoblock;
double pg_sig = PG / (KM_PG + PG);
double couple = 0.35 + 0.65 * (CX43 / CX43_MAX);
double nifblock = IMAX_NIF * CNIF / (IC50_NIF + CNIF);
double cablock = 1.0 - nifblock;
if (cablock < 0.05) cablock = 0.05;
double mgblock = IMAX_MG * pow(CMG, 3.0) / (pow(IC50_MG, 3.0) + pow(CMG, 3.0));
double catarget = (W_TONIC + W_OT * ot_sig + W_PG * pg_sig
+ W_PAR1 * THROMBIN_DRV) * couple * cablock;
dxdt_CAI = KOUT_CAI * (catarget - CAI);
double dcamp = CAMP - CAMPBASE;
if (dcamp < 0.0) dcamp = 0.0;
double camp_brake = CAMP_BRAKE * dcamp / (EC50_CAMP + dcamp);
double ctarget = CONTR_MAX
* pow(CAI, HILL_CONTR) / (pow(EC50_CONTR, HILL_CONTR) + pow(CAI, HILL_CONTR))
* (1.0 - camp_brake) * (1.0 - mgblock);
dxdt_CONTR = KIN_CONTR * (ctarget - CONTR);
// =========================================================================
// 8. Cervical remodelling — the limb tocolysis does not touch
// =========================================================================
double mmpdrive = (0.65 * NFKB + 0.30 * PRW + 0.85 * TERMDRIVE
+ 0.30 * CERV_FRAGILE + 0.20 * (PG / (1.0 + PG)))
* (1.0 - P4_CERV * p4prot);
dxdt_MMP = KIN_MMP * mmpdrive * (MMP_MAX - MMP) - KOUT_MMP * MMP;
dxdt_HA = KIN_HA * (0.90 * TERMDRIVE + 0.45 * PRW + 0.30 * NFKB) * (1.0 - HA)
- KOUT_HA * HA;
// Collagen loss is COOPERATIVE in MMP activity — a low background of MMP
// does not dissolve the cervix, which is why physiologic term shortening is
// gradual while inflamed ripening is abrupt.
double mmpfrac = MMP / MMP_MAX;
dxdt_COLL = -KDEG_COLL * pow(mmpfrac, HILL_MMPCOL) * COLL
+ KREP_COLL * (1.0 - COLL) * (1.0 - 0.8 * NFKB);
// Structural target length: collagen integrity sets the ceiling; hyaluronan
// water imbibition and intrinsic fragility subtract; cerclage/pessary add an
// effective mechanical offset.
double cl_struct = CL_MIN + (CL0 - CL_MIN) * COLL * (1.0 - 0.45 * HA)
* (1.0 - 0.25 * CERV_FRAGILE)
+ CERC_SUPPORT * CERC_EFF + PESS_SUPPORT * PESS_EFF;
// Contractions transiently shorten the cervix (dynamic shortening / funnelling)
double cl_tgt = cl_struct - 2.5 * CONTR / (CONTR_REF + CONTR);
if (cl_tgt < CL_MIN) cl_tgt = CL_MIN;
dxdt_CLEN = KCL * (cl_tgt - CLEN);
// =========================================================================
// 9. Fetal membranes — the other limb tocolysis does not touch
// =========================================================================
double membdeg = KDEG_MEMB * (0.6 * MMP / MMP_MAX + 0.4 * NFKB)
* (1.0 + 0.6 * STRETCH);
dxdt_MEMB = -membdeg * MEMB + KREP_MEMB * (1.0 - MEMB);
double ffndrive = 0.55 * (1.0 - MEMB) + 0.30 * NFKB
+ 0.25 * (1.0 - CLEN / (CL0 + 1e-6));
if (ffndrive < 0.0) ffndrive = 0.0;
dxdt_FFN = KIN_FFN * ffndrive * ffndrive - KOUT_FFN * FFN;
// =========================================================================
// 10. Fetal lung maturation and brain magnesium exposure
// The corticosteroid term is written so the benefit accrues over ~24 h
// and then wanes — the empirical 24 h to 7 d window.
// =========================================================================
double surf_endo = 1.0 / (1.0 + exp(-(GA - GA_SURF50) / GA_SURF_ST));
double acs_eff = EMAX_SURF * BET_F / (EC50_BETF + BET_F);
double excess = (SURF > surf_endo) ? (SURF - surf_endo) : 0.0;
dxdt_SURF = KIN_SURF * (surf_endo * (1.0 + 0.25 * (CORTF - 1.0)) - SURF)
+ KSURF_ACS * acs_eff * (1.0 - SURF)
- ACS_DECAY * excess;
if (SURF > 1.0) dxdt_SURF -= 5.0 * (SURF - 1.0);
dxdt_MGBRAIN = KMG_BRAIN * ((CMG - MG_BASE) - MGBRAIN);
// =========================================================================
// 11. THE DELIVERY HAZARD — a PRODUCT of three limbs
// This multiplicative form is the model's central structural claim: a
// drug acting on one factor cannot push lambda below the product of the
// other two, however completely it silences its own target.
// =========================================================================
double f_contr = CONTR_FLOOR + CONTR / CONTR_REF;
double f_cerv = CERV_FLOOR + CERV_AMP / (1.0 + exp((CLEN - CL_REF) / CL_STEEP));
double f_memb = 1.0 + MEMB_AMP / (1.0 + exp((MEMB - MEMB_REF) / MEMB_STEEP));
double f_term = 1.0 + TERM_AMP * TERMDRIVE;
double ruptured = (MEMB < 0.30 || PPROM_FLAG > 0.5) ? 1.0 : 0.0;
double f_pprom = 1.0 + (PPROM_HAZ - 1.0) * ruptured;
dxdt_CUMHAZ = LAMBDA0 * f_contr * f_cerv * f_memb * f_term * f_pprom;
// Cumulative fetal inflammatory exposure — the currency in which prolonging
// pregnancy inside an infected uterus is PAID for.
dxdt_CUMINFL = IL6;
// =========================================================================
// 12. Drug PK
// =========================================================================
// vaginal micronized progesterone (mg depot -> ng/mL x L plasma amount)
dxdt_P4V_DEP = -KA_P4V * P4V_DEP;
dxdt_P4_C = KA_P4V * P4V_DEP * F_P4V * 1000.0 - (CL_P4 / VD_P4) * P4_C;
// 17-OHPC: absorption-rate-limited (flip-flop) IM depot
dxdt_OHPC_DEP = -KA_OHPC * OHPC_DEP;
dxdt_OHPC_C = KA_OHPC * OHPC_DEP * F_OHPC - (CL_OHPC / VD_OHPC) * OHPC_C;
// atosiban, 2-compartment
dxdt_ATO_C = -(CL_ATO / VC_ATO) * ATO_C - (Q_ATO / VC_ATO) * ATO_C
+ (Q_ATO / VP_ATO) * ATO_P;
dxdt_ATO_P = (Q_ATO / VC_ATO) * ATO_C - (Q_ATO / VP_ATO) * ATO_P;
// nifedipine
dxdt_NIF_GUT = -KA_NIF * NIF_GUT;
dxdt_NIF_C = KA_NIF * NIF_GUT * F_NIF - (CL_NIF / VD_NIF) * NIF_C;
// indomethacin
dxdt_IND_GUT = -KA_IND * IND_GUT;
dxdt_IND_C = KA_IND * IND_GUT * F_IND - (CL_IND / VD_IND) * IND_C;
// magnesium: 2-compartment, RENAL elimination only (CrCl-scaled)
dxdt_MG_C = -(CLMG / VC_MG) * MG_C - (Q_MG / VC_MG) * MG_C + (Q_MG / VP_MG) * MG_P;
dxdt_MG_P = (Q_MG / VC_MG) * MG_C - (Q_MG / VP_MG) * MG_P;
// betamethasone: two IM depots (phosphate ester fast, acetate slow)
dxdt_BET_DP = -KA_BETP * BET_DP;
dxdt_BET_DA = -KA_BETA * BET_DA;
dxdt_BET_C = KA_BETP * BET_DP + KA_BETA * BET_DA - (CL_BET / VD_BET) * BET_C;
dxdt_BET_F = KEQ_BETF * (F_PLAC_BET * CBET - BET_F);
// terbutaline
dxdt_TERB_DEP = -KA_TERB * TERB_DEP;
dxdt_TERB_C = KA_TERB * TERB_DEP - (CL_TERB / VD_TERB) * TERB_C;
[TABLE]
double GAW = GA0_WK + TIME / 7.0;
double CP4o = P4_C / VD_P4;
double CP4LOCo = CP4o * FUP_RATIO;
double COHPCo = 1000.0 * OHPC_C / VD_OHPC;
double CATOo = 1000.0 * ATO_C / VC_ATO;
double CNIFo = 1000.0 * NIF_C / VD_NIF;
double CINDo = 1000.0 * IND_C / VD_IND;
double CMGo = MG_BASE + MG_C / VC_MG;
double CBETo = 1000.0 * BET_C / VD_BET;
double CTERBo = 1000.0 * TERB_C / VD_TERB;
// cumulative probability of having delivered by TIME
double PDELIV = 1.0 - exp(-CUMHAZ);
// biomarker readouts
double IAI = (IL6 > IL6_THRESH) ? 1.0 : 0.0;
double FFNPOS = (FFN > FFN_POS) ? 1.0 : 0.0;
double BISHOPish = 13.0 * (1.0 - CLEN / 40.0);
if (BISHOPish < 0.0) BISHOPish = 0.0;
// safety readouts ---------------------------------------------------------
// magnesium toxicity ladder (mg/dL): 4-7 therapeutic, 9-12 loss of deep
// tendon reflexes, 12-18 respiratory depression, >18 cardiac arrest
double MG_DTR = (CMGo > 9.0) ? 1.0 : 0.0;
double MG_RESP = (CMGo > 12.0) ? 1.0 : 0.0;
// indomethacin ductal constriction — susceptibility is steeply GA-dependent
double DUCTSUS = 1.0 / (1.0 + exp(-(GAW - DUCT_GA50) / DUCT_GA_ST));
double CIND_F = CINDo * F_PLAC_IND;
double DUCTRISK = DUCTSUS * CIND_F / (300.0 + CIND_F);
// betamethasone maternal hyperglycaemia surrogate
double HYPERGLY = CBETo / (4.0 + CBETo);
// nifedipine maternal MAP reduction surrogate
double MAP_DROP = 14.0 * CNIFo / (IC50_NIF + CNIFo);
double B2_LOST = DESENS;
// hazard, exported limb by limb so the decomposition can be inspected -----
double TERMD = 1.0 / (1.0 + exp(-(GAW - GA_RIPEN) / GA_STEEP));
double FCONTR = CONTR_FLOOR + CONTR / CONTR_REF;
double FCERV = CERV_FLOOR + CERV_AMP / (1.0 + exp((CLEN - CL_REF) / CL_STEEP));
double FMEMB = 1.0 + MEMB_AMP / (1.0 + exp((MEMB - MEMB_REF) / MEMB_STEEP));
double RUPT = (MEMB < 0.30 || PPROM_FLAG > 0.5) ? 1.0 : 0.0;
double LAMBDA = LAMBDA0 * FCONTR * FCERV * FMEMB * (1.0 + TERM_AMP * TERMD)
* (1.0 + (PPROM_HAZ - 1.0) * RUPT);
[CAPTURE] @annotated
GAW : weeks, gestational age
CP4o : ng/mL, plasma progesterone (exogenous component)
CP4LOCo : ng/mL-eq, LOCAL endometrial progesterone (first uterine pass)
COHPCo : ng/mL, 17-OHPC plasma concentration
CATOo : ng/mL, atosiban plasma concentration
CNIFo : ng/mL, nifedipine plasma concentration
CINDo : ng/mL, indomethacin plasma concentration
CMGo : mg/dL, serum magnesium
CBETo : ng/mL, maternal betamethasone
CTERBo : ng/mL, terbutaline
PDELIV : -, cumulative probability of delivery by this time
LAMBDA : 1/day, instantaneous delivery hazard
FCONTR : -, contraction limb of the hazard
FCERV : -, cervical limb of the hazard
FMEMB : -, membrane limb of the hazard
TERMD : -, physiologic term drive
IAI : 0/1, intra-amniotic inflammation (IL-6 > 2.6 ng/mL)
FFNPOS : 0/1, fetal fibronectin positive (> 50 ng/mL)
BISHOPish: -, Bishop-like cervical score derived from CL
MG_DTR : 0/1, magnesium above the DTR-loss threshold
MG_RESP : 0/1, magnesium above the respiratory-depression threshold
DUCTSUS : -, gestational-age-dependent ductal susceptibility
DUCTRISK : -, indomethacin ductal-constriction risk index
HYPERGLY : -, maternal hyperglycaemia index (betamethasone)
MAP_DROP : mmHg, nifedipine-associated MAP reduction
B2_LOST : -, fraction of beta2 responsiveness lost (tachyphylaxis)
[ENV]
## ============================================================================
## R DRIVER CODE ($ENV — evaluated when the model is loaded with mread())
## ----------------------------------------------------------------------------
## Usage:
## library(mrgsolve); library(dplyr)
## mod <- mread("ptl_mrgsolve_model.R")
## # bring the helpers below into the global environment:
## list2env(as.list(mrgsolve::env_get(mod)), .GlobalEnv)
##
## out <- PTL_simulate_scenarios(mod)
## PTL_ledger(mod) # days gained vs OUTCOME gained
## PTL_limb_decomposition(mod) # why tocolysis is arithmetically bounded
## PTL_acs_window(mod) # the 24 h - 7 d steroid window
## PTL_infection_paradox(mod) # days up, outcome down
## PTL_tocolytic_head_to_head(mod)
## ============================================================================
## ---------------------------------------------------------------------------
## DOSING REGIMEN BUILDERS
## ---------------------------------------------------------------------------
# --- Prevention -------------------------------------------------------------
#' Vaginal micronized progesterone 200 mg capsule nightly
ev_vaginal_P4 <- function(start_d = 0, stop_d = 126, dose_mg = 200) {
mrgsolve::ev(time = start_d, amt = dose_mg, cmt = "P4V_DEP",
ii = 1, addl = max(0, floor(stop_d - start_d)))
}
#' 17-OHPC 250 mg IM weekly (classically 16+0 to 36+0 weeks)
ev_17OHPC <- function(start_d = 0, stop_d = 112, dose_mg = 250) {
n <- max(0, floor((stop_d - start_d) / 7))
mrgsolve::ev(time = start_d, amt = dose_mg, cmt = "OHPC_DEP", ii = 7, addl = n)
}
# --- Acute tocolysis --------------------------------------------------------
#' Atosiban: 6.75 mg IV bolus, then 300 ug/min x 3 h, then 100 ug/min x 45 h
#' (bolus + two zero-order infusions)
ev_atosiban <- function(start_d = 0) {
bolus <- mrgsolve::ev(time = start_d, amt = 6.75, cmt = "ATO_C")
# 300 ug/min = 432 mg/day for 3 h (0.125 d) -> 54 mg
hi <- mrgsolve::ev(time = start_d, amt = 54, cmt = "ATO_C", rate = 432)
# 100 ug/min = 144 mg/day for 45 h (1.875 d) -> 270 mg
lo <- mrgsolve::ev(time = start_d + 0.125, amt = 270, cmt = "ATO_C", rate = 144)
c(bolus, hi, lo)
}
#' Nifedipine IR: 20 mg PO loading, then 20 mg q6h
ev_nifedipine <- function(start_d = 0, hours = 48) {
n <- max(0, floor(hours / 6) - 1)
mrgsolve::ev(time = start_d, amt = 20, cmt = "NIF_GUT", ii = 0.25, addl = n)
}
#' Indomethacin: 100 mg PR loading, then 25 mg q6h for <= 48 h
ev_indomethacin <- function(start_d = 0, hours = 48) {
load <- mrgsolve::ev(time = start_d, amt = 100, cmt = "IND_GUT")
n <- max(0, floor(hours / 6) - 1)
main <- mrgsolve::ev(time = start_d + 0.25, amt = 25, cmt = "IND_GUT",
ii = 0.25, addl = n)
c(load, main)
}
#' MgSO4: `load_g` grams IV over 30 min, then `rate_g_h` g/h.
#' Model units: 1 g MgSO4 -> 10 (mg/dL x L) in the central compartment, so a
#' 4 g load into VC_MG = 16 L raises serum Mg by ~2.5 mg/dL.
ev_magnesium <- function(start_d = 0, hours = 24, load_g = 4, rate_g_h = 1) {
u <- 10
load <- mrgsolve::ev(time = start_d, amt = load_g * u, cmt = "MG_C",
rate = load_g * u / (0.5 / 24))
maint <- mrgsolve::ev(time = start_d + 0.5 / 24,
amt = rate_g_h * u * hours, cmt = "MG_C",
rate = rate_g_h * u * 24)
c(load, maint)
}
#' Betamethasone 12 mg IM x 2, 24 h apart (phosphate + acetate, 1:1)
ev_betamethasone <- function(start_d = 0, n_doses = 2, dose_mg = 12) {
p <- mrgsolve::ev(time = start_d, amt = dose_mg / 2, cmt = "BET_DP",
ii = 1, addl = n_doses - 1)
a <- mrgsolve::ev(time = start_d, amt = dose_mg / 2, cmt = "BET_DA",
ii = 1, addl = n_doses - 1)
c(p, a)
}
#' Terbutaline 0.25 mg SC q4h (historical; FDA boxed warning against
#' prolonged use for tocolysis)
ev_terbutaline <- function(start_d = 0, hours = 24) {
n <- max(0, floor(hours / 4) - 1)
mrgsolve::ev(time = start_d, amt = 0.25, cmt = "TERB_DEP",
ii = 4 / 24, addl = n)
}
#' A null event (no drug) — needed because mrgsim_e() wants an event object
ev_none <- function() mrgsolve::ev(time = 0, amt = 0, cmt = "P4V_DEP")
## ---------------------------------------------------------------------------
## GESTATIONAL-AGE-DEPENDENT NEONATAL OUTCOME MODEL
## ---------------------------------------------------------------------------
#' Neonatal outcome probabilities given gestational age at delivery, the
#' lung-maturity index, the fetal inflammatory state and magnesium exposure.
PTL_neonatal_outcomes <- function(ga_wk, surf = 0.2, firs = 0, mgbrain = 0,
transfer = 0, coamox = 0, p = list()) {
d <- list(RDS_GA50 = 30.2, RDS_STEEP = 1.9, RDS_SURF_B = 2.0,
IVH_GA50 = 26.5, IVH_STEEP = 1.5,
DEATH_GA50 = 24.3, DEATH_STEEP = 1.4,
NEC_GA50 = 26.0, NEC_STEEP = 2.0,
CP_BASE = 0.020, EMAX_MGN = 0.45, EC50_MGN = 1.2,
FIRS_ON_CP = 2.00, FIRS_ON_SEP = 0.45, SEP_BASE = 0.020,
FIRS_ON_IVH = 1.00, FIRS_ON_NEC = 1.40, FIRS_ON_DEATH = 0.60)
d[names(p)] <- p
inv <- function(x) 1 / (1 + exp(-x))
# A fetal inflammatory response is not a spectator: it raises IVH, NEC,
# sepsis, death and cerebral palsy. This is the term that makes buying days
# inside an infected amniotic cavity a LOSING trade.
rds <- inv(-(ga_wk - d$RDS_GA50) / d$RDS_STEEP - d$RDS_SURF_B * surf)
ivh <- inv(-(ga_wk - d$IVH_GA50) / d$IVH_STEEP) * (1 - 0.45 * surf) *
(1 + d$FIRS_ON_IVH * firs)
nec <- inv(-(ga_wk - d$NEC_GA50) / d$NEC_STEEP) * (1 + d$FIRS_ON_NEC * firs) *
(1 + 0.35 * coamox)
death <- inv(-(ga_wk - d$DEATH_GA50) / d$DEATH_STEEP) * (1 - 0.22 * surf) *
(1 + d$FIRS_ON_DEATH * firs) * (1 - 0.20 * transfer)
sepsis <- d$SEP_BASE + d$FIRS_ON_SEP * firs
cp_ga <- d$CP_BASE * (1 + 2.4 * inv(-(ga_wk - 28.5) / 2.0))
mg_eff <- d$EMAX_MGN * mgbrain / (d$EC50_MGN + mgbrain)
cp <- cp_ga * (1 + d$FIRS_ON_CP * firs) * (1 - mg_eff)
comp <- 1 - (1 - rds) * (1 - ivh) * (1 - nec) * (1 - sepsis)
data.frame(GA_wk = ga_wk, RDS = rds, severe_IVH = ivh, NEC = nec,
sepsis = sepsis, death = death, CP = cp, composite = comp)
}
## ---------------------------------------------------------------------------
## SCENARIOS
## ---------------------------------------------------------------------------
#' Eleven prebuilt scenarios spanning prevention, acute tocolysis, PPROM and
#' the infection case in which tocolysis is actively harmful.
PTL_scenarios <- function() {