Repository navigation
Expand file tree
/
Copy pathlipo_mrgsolve_model.R
More file actions
1285 lines (1166 loc) · 61.5 KB
/
Copy pathlipo_mrgsolve_model.R
File metadata and controls
1285 lines (1166 loc) · 61.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
## =============================================================================
## lipo_mrgsolve_model.R
## Lipodystrophy syndromes (CGL · AGL · FPLD · APL) — Quantitative Systems
## Pharmacology model for mrgsolve
##
## 53 ODE compartments. Time unit = DAYS (a 4-hour metreleptin half-life and
## a 10-year fibrosis trajectory have to live in the same system).
##
## ---------------------------------------------------------------------------
## THE FOUR STRUCTURAL COMMITMENTS
## ---------------------------------------------------------------------------
##
## (1) ONE LESION, TWO DEFICITS, AND THE DRUGS ARE SORTED BY WHICH ONE THEY
## TOUCH.
##
## Adipose tissue is simultaneously a MECHANICAL BUFFER (the only tissue
## that can esterify a large triglyceride flux cheaply) and the ENDOCRINE
## REPORTER of how large that buffer is. Losing adipocytes therefore
## produces two deficits that are mechanistically independent:
##
## capacity deficit : J_ov = max(0, J_in - Sum(J_st,i) - OX)
## signal deficit : DEF = 1 - SD(L_effective)
##
## and the signal deficit RAISES J_in through hyperphagia, so the two
## deficits multiply rather than add. Every drug in this file is
## classified by which term it enters:
##
## metreleptin -> SD (signal) ... and hence J_in
## pioglitazone -> C_cap (capacity), gated by preadipocyte substrate
## APOC3 / ANGPTL3 / fibrate / omega-3 -> plasma TG clearance only
## insulin / metformin / SGLT2i -> glucose only
## very-low-fat diet -> J_in directly
##
## (2) THE GL/PL DISSOCIATION IS NOT WRITTEN DOWN. IT IS DERIVED.
##
## There is no parameter in this model called "metreleptin works better in
## generalised lipodystrophy". There is one saturating transducer
##
## S(L) = L^h / (L^h + EC50^h), EC50 = 4 ng/mL, h = 2
##
## applied to total effective leptin. Because S saturates, the SAME dose
## produces the same receptor occupancy in every patient, while the
## achievable change in S is bounded by 1 - S(L_baseline) -- a property of
## the PATIENT, not of the dose. Occupancy is captured as OCC and the
## deficit as DEF so that the product structure can be read off the
## output. Scenarios 7 / 8 / 9 exist to make the claim falsifiable:
## identical mg/kg in two phenotypes, then a dose escalation that raises
## OCC and should NOT rescue the effect if the model is right.
##
## (3) PLASMA TRIGLYCERIDE IS TIER 5 OF A CASCADE, NOT THE DISEASE.
##
## Overflow is partitioned to liver, muscle, plasma, pancreas, heart and
## kidney with fixed fractions. Plasma TG is therefore a MARKER of J_ov,
## and an agent that removes it (APOC3 ASO) must leave hepatic fat and
## HbA1c nearly unchanged. That decoupling is a prediction, not an
## assumption -- scenario 15 is the test, and if hepatic fat fell with
## plasma TG the tier structure would be wrong.
##
## (4) ONLY FIBROSIS IS IRREVERSIBLE.
##
## Every other state in this model relaxes back when the driver is
## removed (scenario 10 withdraws metreleptin and everything returns).
## FIB has a regression rate 17x slower than its formation rate, so early
## versus late initiation (scenario 23) changes the CEILING of recovery
## rather than its rate. No bistability is claimed anywhere: the model
## has no positive feedback loop with gain > 1, and the withdrawal
## scenario is included precisely so that claim can be checked.
##
## ---------------------------------------------------------------------------
## WHAT THIS MODEL IS CALIBRATED AGAINST
## ---------------------------------------------------------------------------
## Baseline phenotypes (the model is run to its own steady state for each
## phenotype -- see baseline() -- rather than having initial conditions typed
## in by hand):
##
## congenital generalised (CGL): leptin 0.4-1.5 ng/mL, HFF 20-35%,
## TG 800-3000 mg/dL, HbA1c 9-11%, insulin 40-100 uU/mL, ALT 60-120 U/L
## (Garg 2004 NEJM; Agarwal 2003; Simha 2003 JCEM; Akinci 2019)
## familial partial (FPLD2, Dunnigan): leptin 4-12 ng/mL, HFF 10-20%,
## TG 250-600 mg/dL, HbA1c 7.5-8.5% (Vigouroux 2011; Guillin-Amarelle 2016)
## healthy female 60 kg: leptin 10-15 ng/mL, HFF ~3%, TG ~90 mg/dL,
## HbA1c ~5.1%, insulin ~8 uU/mL
##
## Treatment effects the model is asked to reproduce (see summarise_all()):
## metreleptin in GL HbA1c -1.7 to -2.5%, TG -40 to -60%,
## hepatic fat -30 to -50% relative, insulin dose
## requirement often halved
## (Oral 2002 NEJM; Petersen 2002 JCI; Javor 2005;
## Chan 2011; Diker-Cohen 2015 JCEM;
## Brown 2018 Endocrine Reviews / MEASuRE cohort)
## metreleptin in PL HbA1c ~ -0.6% overall, benefit concentrated in
## the low-leptin stratum
## (Diker-Cohen 2015; Oral 2019 NEJM combined analysis)
## volanesorsen in FPLD TG -70 to -88%, hepatic fat change small
## (BROADEN, Oral 2019/Prohaska 2023)
## pioglitazone in FPLD modest fat redistribution + glycaemic gain;
## no adipogenic gain where no preadipocytes exist
##
## Parameters that are honest guesses rather than fitted values are marked
## "(struct)" in the annotation. They set the SHAPE of a relationship whose
## existence is documented but whose coefficient is not: e.g. the fraction of
## overflow going to pancreas, or the slope of proteinuria on insulin.
##
## ---------------------------------------------------------------------------
## UNITS
## ---------------------------------------------------------------------------
## time days
## metreleptin micrograms (so amount/L = ng/mL, the clinical unit)
## small molecules mg, L
## adipose depots kg of triglyceride
## liver / muscle TG grams
## plasma TG, glucose mg/dL
## insulin microunits/mL
## leptin ng/mL
## normalised states 1 = the healthy steady state (IMCLN, PLIP, MYOL, KLIP)
##
## Author: QSP Disease Model Library (Claude Code Routine)
## Licence: see repository LICENSE. EDUCATIONAL / RESEARCH USE ONLY —
## not validated for clinical or regulatory use.
## =============================================================================
library(mrgsolve)
lipo_code <- r"---(
$PROB
# Lipodystrophy syndromes QSP model (53 ODEs)
$PARAM @annotated
// ================================================================= PHENOTYPE
// The phenotype is a CAPACITY DISTRIBUTION plus a secretion coefficient.
// There is deliberately no categorical "disease type" switch.
FC1 : 1.00 : femorogluteal SC capacity, fraction of normal
FC2 : 1.00 : upper-trunk/facial SC capacity, fraction of normal
FCV : 1.00 : visceral capacity, fraction of normal
PREAD : 1.00 : preadipocyte/expandability reserve, 0-1 (gates PPARg effect)
LEPKM : 1.00 : multiplier on leptin secretion per kg fat (1 = normal biology)
LEAN0 : 45 : lean body mass (kg)
BASLIP : 0 : excess basal lipolysis (PLIN1/CIDEC phenotypes), 0-1 (struct)
// ============================================================ CAPACITY (kg)
CAP1N : 12.0 : normal femorogluteal TG ceiling (kg)
CAP2N : 5.0 : normal upper-trunk/facial TG ceiling (kg)
CAPVN : 3.0 : normal visceral TG ceiling (kg)
FATN : 14.0 : reference healthy total fat mass (kg), scales LPL and leptin
// ================================================== ENERGY INTAKE AND INFLUX
EIBASE : 40.0 : weight-maintenance intake per kg lean mass (kcal/kg/d)
FFAT : 0.35 : habitual fraction of energy as fat
FCHO : 0.45 : habitual fraction of energy as carbohydrate
HMAX : 0.45 : maximal fractional hyperphagia at full signal deficit
TAUEI : 3.0 : intake adaptation time constant (d)
EICLAMP : 0 : >0 fixes intake (kcal/d) = pair-feeding clamp
DIETT : 1e6 : day a very-low-fat diet starts
DIETFAT : 0.15 : fat fraction of energy on that diet (CHO takes the balance)
KDNL : 0.030 : de novo lipogenesis yield per g carbohydrate at full insulin drive
KINS50 : 20.0 : insulin giving half-maximal lipogenic drive (uU/mL)
KDNLSEL : 0.60 : extra DNL drive from selective hepatic insulin resistance (struct)
KLIPSC : 0.0100 : subcutaneous adipose TG turnover (1/d) -- the futile cycle
KLIPVI : 0.0120 : visceral adipose TG turnover (1/d)
// ============================================== BUFFER (ESTERIFICATION) TERMS
KEST : 15.0 : esterification Vmax per kg of depot capacity (g/d/kg)
HRFLOOR : 0.10 : irreducible turnover headroom of a full depot (struct)
KMJ : 60.0 : influx giving half-maximal storage flux (g/d)
KOXW : 2.00 : EXTRAhepatic fat oxidation coefficient (g/d per kg^0.75)
KOXAD : 0.35 : adaptive rise in oxidation with ectopic load (struct)
TAUOXAD : 20.0 : time constant of that adaptation (d)
// ====================================================== OVERFLOW PARTITIONING
FLIV : 0.35 : baseline fraction of overflow to liver
KFLIVA : 0.30 : extra hepatic share as adipose disappears (competitive sinks)
FMUS : 0.15 : fraction of overflow to skeletal muscle
FPL : 0.15 : fraction of overflow appearing directly in plasma
FPANC : 0.03 : fraction of overflow to pancreatic islets (struct)
FMYO : 0.04 : fraction of overflow to myocardium (struct)
FKID : 0.03 : fraction of overflow to kidney (struct)
// ================================================================ LIVER LIPID
LIVLEAN : 1400 : lean liver mass (g), sets hepatic fat fraction
FHNEFA : 0.20 : splanchnic share of the systemic NEFA flux taken by the liver
FHVIS : 0.55 : PORTAL share of visceral lipolysis delivered to the liver
SCFATN : 11.6 : healthy SUBCUTANEOUS fat mass (kg) -- the liver's only competitor
VMVLDL : 45.0 : maximal VLDL-TG secretion (g/d)
KMLTG : 68.0 : hepatic TG giving half-maximal VLDL secretion (g)
VMOX : 19.0 : maximal hepatic FA oxidation (g/d) -- extrahepatic oxidation
KMOX : 68.0 : hepatic TG for half-maximal oxidation (g) -- is in KOXW
KLVLDL : 0.00 : direct leptin suppression of VLDL secretion -- SET TO ZERO, see \$ODE
KLOX : 0.25 : hepatic FA-oxidation LOSS per unit leptin signal lost
KLDNL : 0.20 : de novo lipogenesis gain per unit leptin signal lost
KOM3VLDL : 0.18 : omega-3 suppression of VLDL-TG secretion
KEVVLDL : 0.05 : ANGPTL3-blockade suppression of VLDL-TG secretion
KFIBOX : 0.50 : PPARa (fibrate) gain on hepatic FA oxidation
// =============================================================== PLASMA LIPID
VPDL : 30.0 : plasma volume (dL)
VMTG : 1730 : maximal LPL-mediated TG clearance (mg/dL/d)
KMTGP : 150 : plasma TG at half-maximal LPL clearance (mg/dL)
KLINTG : 0.90 : non-saturable TG clearance (1/d)
LPLFLOOR : 0.45 : LPL activity retained with no adipose tissue (muscle, heart)
KAP3 : 2.50 : LPL gain per unit APOC3 knockdown
KANG : 1.20 : LPL gain per unit ANGPTL3 knockdown
KFIBL : 0.50 : LPL gain from PPARa activation
TGTHR : 1000 : pancreatitis threshold for time-above-threshold (mg/dL)
KPANC : 0.0015 : maximal pancreatitis hazard (1/d) (struct)
// ================================================================ NEFA / IMCL
NEFA0 : 0.45 : reference plasma NEFA (mmol/L)
KLIPB : 1.20 : NEFA gain from excess basal lipolysis
KINSN : 0.60 : insulin suppression of NEFA (per baseline-insulin unit)
TAUNEFA : 0.50 : NEFA turnover time constant (d)
IMCLN0 : 39.5 : reference intramyocellular lipid pool (g)
KIMOX : 0.15 : IMCL disposal (1/d)
TAUECT : 30.0 : time constant of pancreas/heart/kidney lipid tracking (d)
JOVREF : 39.5 : healthy reference overflow flux (g/d), normalises ectopic pools
KINSACT : 8.00 : insulin action (level x sensitivity) at half-maximal esterification
// ==================================================== INSULIN RESISTANCE INDEX
AHFF : 0.220 : IR slope on hepatic fat fraction (per %)
BIM : 0.150 : IR slope on normalised IMCL
CINFL : 0.400 : IR slope on hepatic inflammation (struct)
CLEP : 2.000 : IR reduction per unit leptin signal (direct peripheral action)
DADP : 1.000 : IR reduction per unit relative adiponectin
SDREF : 0.900 : leptin signal of the healthy reference subject
ADPN0 : 7.50 : healthy adiponectin (ug/mL)
ADPNFLR : 0.50 : non-adipose adiponectin floor (ug/mL)
IRMIN : 0.50 : floor on the IR index
TAUIR : 7.00 : lag of tissue insulin resistance behind lipid load (d)
// ============================================================ GLUCOSE-INSULIN
VGDL : 96.0 : glucose distribution volume (dL)
HGP0 : 2400 : reference hepatic glucose output (mg/dL/d)
KHIR : 0.38 : HGP gain per unit excess IR
KHINS : 0.080 : HGP suppression per uU/mL insulin
KGZ : 4.00 : insulin-independent glucose disposal (1/d)
KGI : 1.53 : insulin-dependent disposal (1/d per uU/mL)
KREN : 8.00 : renal glucose excretion above threshold (1/d)
GRENTHR : 180 : renal glucose threshold (mg/dL)
SMAX : 960 : maximal beta-cell secretory rate (uU/mL/d)
G50 : 110 : glucose for half-maximal secretion (mg/dL)
KCOMP : 0.35 : compensatory secretory gain per unit excess IR
KINSCL : 48.0 : insulin elimination (1/d) -- effective, smoothed
KIEXT : 0.15 : loss of hepatic insulin extraction per unit excess IR
KBREC : 0.0020 : beta-cell functional recovery (1/d)
KBLIP : 0.0010 : beta-cell loss per unit excess islet lipid (1/d)
KBGLU : 0.0015 : beta-cell loss per 100 mg/dL glucose above 180 (1/d)
GTOXTHR : 180 : glucose above which glucotoxicity operates (mg/dL)
TAUA1C : 35.0 : HbA1c integration time constant (d)
MGFAC : 1.10 : mean-glucose to model-glucose ratio
// ================================================ LIVER INJURY AND FIBROSIS
ALT0 : 18.0 : healthy ALT (U/L)
KALT : 2.60 : ALT rise per % hepatic fat
TAUALT : 10.0 : ALT turnover (d)
K50INFL : 15.0 : hepatic fat fraction for half-maximal inflammation (%)
TAUINFL : 30.0 : inflammation turnover (d)
KFIB : 0.0100 : fibrogenesis rate above the inflammation threshold (1/d)
INFLTHR : 0.35 : inflammation threshold for net fibrogenesis
KFREG : 0.00015 : fibrosis regression (1/d) -- 67x slower than formation
// ==================================================== KIDNEY / ANDROGEN / QoL
KPR : 0.15 : proteinuria slope on excess IR (g/d)
KPRG : 0.20 : proteinuria slope per 100 mg/dL glucose above 140 (g/d)
KPRL : 0.25 : proteinuria slope on renal lipid (g/d) (struct)
TAUPROT : 60.0 : proteinuria adaptation (d)
KEGFR : 0.005 : eGFR loss per g/d proteinuria (mL/min/1.73/d)
KAND : 0.55 : androgen index gain per relative insulin excess (struct)
TAUAND : 30.0 : androgen index adaptation (d)
// ==================================================== LEPTIN AXIS / TRANSDUCER
LEPK : 0.60 : leptin secretion per kg effective fat (ng/mL/kg)
HYPMAX : 0.60 : extra leptin per kg from adipocyte hypertrophy at full fill
LEPFLOOR : 0.15 : non-adipose leptin floor (ng/mL)
TAULEP : 0.50 : endogenous leptin turnover (d)
LEC50 : 4.00 : leptin for half-maximal hypothalamic signal (ng/mL)
LHILL : 2.00 : Hill coefficient of the leptin transducer
KSOCS : 3.00 : SOCS3/PTP1B brake gain on signal sustained ABOVE the healthy level
SREFB : 0.88 : healthy hypothalamic signal above which the brake engages
TAUSOCS : 14.0 : brake adaptation time constant (d)
TAUSLAG : 2.00 : lag of peripheral leptin actions behind plasma signal (d)
KDLEP : 6.00 : apparent Kd for the receptor-occupancy readout (ng/mL)
// ====================================================== METRELEPTIN PK AND ADA
KAMET : 8.00 : metreleptin SC absorption (1/d)
FMET : 0.80 : metreleptin SC bioavailability
VMETC : 35.0 : metreleptin central volume (L)
CLMET : 145 : metreleptin clearance (L/d) -> t1/2 ~ 4 h
QMET : 12.0 : metreleptin distribution clearance (L/d)
VMETP : 20.0 : metreleptin peripheral volume (L)
NEUT : 0 : neutralising-ADA switch (0/1) -- binding-only ADA has no PD effect
KNABON : 0.030 : neutralising titre formation (1/d)
KNABOFF : 0.0020 : neutralising titre decay (1/d)
NABMAX : 20.0 : maximal neutralising capacity (dimensionless)
CNAB50 : 10.0 : metreleptin concentration driving titre formation (ng/mL)
// ================================================================ PIOGLITAZONE
KAPIO : 6.00 : pioglitazone absorption (1/d)
VPIO : 15.0 : pioglitazone volume (L)
CLPIO : 20.8 : pioglitazone clearance (L/d), active-metabolite adjusted
EC50PIO : 700 : pioglitazone concentration for half-maximal PPARg effect (ng/mL)
KCAPM : 0.45 : maximal fractional capacity gain at full PPARg x full substrate
TAUCAP : 45.0 : adipogenic capacity expansion time constant (d)
KADPPIO : 6.00 : adiponectin gain from PPARg activation (ug/mL)
TAUADP : 20.0 : adiponectin turnover (d)
// ============================================== VOLANESORSEN (APOC3 ASO) PK/PD
KAVOL : 1.00 : volanesorsen SC absorption (1/d)
VVOLC : 10.0 : volanesorsen central volume (L)
CLVOL : 2.00 : volanesorsen plasma clearance (L/d)
QVOL : 1.50 : plasma-to-tissue clearance (L/d)
VVOLT : 5.00 : tissue (hepatic) volume (L)
KOUTVT : 0.030 : tissue elimination (1/d)
KOUT3 : 0.050 : APOC3 turnover (1/d) -> ~14 d onset
EMAX3 : 0.85 : maximal APOC3 knockdown
EC503 : 1.50 : tissue concentration for half-maximal knockdown (mg/L)
KFIB3 : 0.35 : fibrate-mediated reduction of APOC3 production
// ==================================================== FENOFIBRATE / OMEGA-3
KAFIB : 4.00 : fenofibrate absorption (1/d)
VFIB : 40.0 : fenofibrate volume (L)
CLFIB : 20.0 : fenofibrate clearance (L/d)
EC50FIB : 5.00 : fibrate concentration for half-maximal PPARa effect (mg/L)
OM3DOSE : 0 : omega-3 4 g/d switch (0/1)
OM3T : 1e6 : day omega-3 starts
KOM3OX : 0.35 : omega-3 gain on hepatic FA oxidation (its main action)
TAUOM3 : 14.0 : omega-3 tissue enrichment time constant (d)
// ====================================================== EVINACUMAB (ANGPTL3)
VEVI : 5.00 : evinacumab volume (L)
CLEVI : 0.35 : evinacumab clearance (L/d)
KOUTA : 0.10 : ANGPTL3 turnover (1/d)
EMAXA : 0.90 : maximal ANGPTL3 suppression
EC50A : 10.0 : evinacumab concentration for half-maximal suppression (mg/L)
// ====================================================== GLUCOSE-DIRECTED DRUGS
KAMET2 : 8.00 : metformin absorption (1/d)
VMETF : 600 : metformin volume (L)
CLMETF : 600 : metformin clearance (L/d)
EMAXMETF : 0.25 : maximal metformin suppression of HGP
EC50METF : 1.50 : metformin concentration for half-maximal effect (mg/L)
INSDOSE : 0 : exogenous insulin (U/d) as a continuous requirement
INST : 1e6 : day exogenous insulin starts
KAI : 4.00 : insulin depot absorption (1/d)
SCALEI : 24.0 : uU/mL per U/d of exogenous insulin at steady state
KAGLP : 0.50 : GLP-1 RA absorption (1/d)
VGLP : 12.0 : GLP-1 RA volume (L)
CLGLP : 1.19 : GLP-1 RA clearance (L/d) -> t1/2 ~ 7 d
EMAXGLP : 0.20 : maximal fractional intake reduction by GLP-1 RA
EC50GLP : 100 : GLP-1 RA concentration for half-maximal effect (ng/mL)
KASGL : 8.00 : SGLT2i absorption (1/d)
VSGL : 120 : SGLT2i volume (L)
CLSGL : 20.0 : SGLT2i clearance (L/d)
UGEMAX : 729 : maximal glucosuria (mg/dL/d)
EC50SGL : 0.15 : SGLT2i concentration for half-maximal glucosuria (mg/L)
GSGLTHR : 100 : renal glucose threshold for SGLT2i-driven excretion (mg/dL)
$INIT @annotated
// ---- metreleptin ----------------------------------------------------------
MSC : 0 : metreleptin SC depot (ug)
MCEN : 0 : metreleptin central (ug)
MPER : 0 : metreleptin peripheral (ug)
NAB : 0 : neutralising anti-metreleptin capacity (dimensionless)
// ---- other drugs ----------------------------------------------------------
PGUT : 0 : pioglitazone gut (mg)
PCEN : 0 : pioglitazone central (mg)
VSC : 0 : volanesorsen SC depot (mg)
VCEN : 0 : volanesorsen central (mg)
VTIS : 0 : volanesorsen tissue (mg)
APOC3 : 1 : APOC3 level (fraction of normal)
FGUT : 0 : fenofibrate gut (mg)
FCEN : 0 : fenofibrate central (mg)
EVC : 0 : evinacumab central (mg)
ANG3 : 1 : ANGPTL3 level (fraction of normal)
OM3 : 0 : omega-3 tissue enrichment index (0-1)
MTGUT : 0 : metformin gut (mg)
MTCEN : 0 : metformin central (mg)
ISCD : 0 : exogenous insulin depot (U)
INSX : 0 : exogenous insulin concentration (uU/mL)
GSC : 0 : GLP-1 RA SC depot (mg)
GCEN : 0 : GLP-1 RA central (mg)
SGGUT : 0 : SGLT2i gut (mg)
SGCEN : 0 : SGLT2i central (mg)
// ---- leptin axis and intake ----------------------------------------------
LEPEN : 12.0 : endogenous plasma leptin (ng/mL)
SOCS : 0.5 : hypothalamic desensitisation state
SLAG : 0.9 : lagged leptin signal driving peripheral actions
EIS : 1400 : energy intake state (kcal/d)
// ---- adipose depots and capacity -----------------------------------------
ASC1 : 8.0 : femorogluteal SC TG (kg)
ASC2 : 3.5 : upper-trunk/facial SC TG (kg)
AVIS : 2.0 : visceral TG (kg)
CAPM : 1.0 : capacity multiplier (PPARg-driven adipogenesis)
OXAD : 0 : adaptive oxidation state
// ---- ectopic lipid --------------------------------------------------------
LTG : 45.0 : liver TG (g)
IMCL : 12.0 : intramyocellular lipid (g)
PLIP : 1.0 : islet lipid (normalised)
MYOL : 1.0 : myocardial lipid (normalised)
KLIP : 1.0 : renal lipid (normalised)
// ---- plasma lipid ---------------------------------------------------------
PTG : 90.0 : plasma triglyceride (mg/dL)
NEFA : 0.45 : plasma NEFA (mmol/L)
ADPN : 8.0 : adiponectin (ug/mL)
// ---- glucose axis ---------------------------------------------------------
IRX : 1.0 : lagged insulin resistance index
GLU : 90.0 : plasma glucose (mg/dL)
INS : 8.0 : endogenous insulin (uU/mL)
BCF : 1.0 : beta-cell functional mass (fraction)
A1C : 5.1 : HbA1c (%)
// ---- organ injury ---------------------------------------------------------
ALT : 18.0 : ALT (U/L)
INFL : 0.17 : hepatic inflammation (normalised)
FIB : 0 : fibrosis stage (0-4)
PROT : 0.05 : proteinuria (g/d)
EGFR : 105 : eGFR (mL/min/1.73 m2)
ANDX : 1.0 : androgen index (1 = normal)
TAT : 0 : cumulative time with TG above threshold (d)
PANCH : 0 : cumulative pancreatitis hazard
$GLOBAL
#define HFF (100.0*LTG/(LTG + LIVLEAN))
#define FATM (ASC1 + ASC2 + AVIS)
#define EFFAT (ASC1 + ASC2 + 0.30*AVIS)
// Shared scalars set in $MAIN (constant within a record)
double C1, C2, CV, EI0, DIETFATN, DIETCHON;
// Leptin transducer: saturating Hill on total effective leptin, then the
// SOCS3/PTP1B brake. This single function is the only place where leptin
// becomes a signal, and it is deliberately the same function for endogenous
// and exogenous leptin.
double lep_signal(double L, double ec50, double h, double socs,
double ksocs, double srefb) {
double Lp = pow(fmax(L, 1e-8), h);
double S = Lp/(Lp + pow(ec50, h));
// The SOCS3/PTP1B brake is a DEVIATION mechanism: it engages only for signal
// held above the healthy level, which is why it desensitises obesity without
// taxing a normal subject -- and why it cannot be invoked to explain the
// lipodystrophy phenotype, where the signal is far below srefb.
return S/(1.0 + ksocs*fmax(0.0, socs - srefb));
}
$MAIN
// Capacity ceilings. CAPM (the PPARg state) multiplies only the SUBCUTANEOUS
// depots: thiazolidinediones expand subcutaneous adipose, they do not create
// visceral capacity.
C1 = CAP1N*FC1*CAPM;
C2 = CAP2N*FC2*CAPM;
CV = CAPVN*FCV;
// Weight-maintenance intake scales with lean mass, not with total weight:
// lipodystrophy patients have little fat but normal-to-high lean mass.
EI0 = EIBASE*LEAN0*exp(ETA(1));
DIETFATN = DIETFAT;
DIETCHON = FFAT + FCHO - DIETFAT;
$ODE
// ===========================================================================
// 1. LEPTIN: SECRETION, EXOGENOUS DRUG, TRANSDUCER, DEFICIT
// ===========================================================================
double fill1 = ASC1/(C1 + 1e-9);
double fill2 = ASC2/(C2 + 1e-9);
double fillv = AVIS/(CV + 1e-9);
double fillm = fmin(1.0, (ASC1 + ASC2)/(C1 + C2 + 1e-9));
// Leptin per kg rises with adipocyte size (fractional fill), which is why the
// few remaining, maximally distended adipocytes of a CGL patient secrete more
// per gram than a healthy adipocyte does.
double lep_ss = LEPKM*LEPK*EFFAT*(1.0 + HYPMAX*fillm) + LEPFLOOR;
double CMET = MCEN/VMETC; // total metreleptin (ng/mL)
double FU = 1.0/(1.0 + NAB); // free fraction left by neutralising ADA
double LTOT = LEPEN + CMET; // what a clinical assay measures
double LEFF = LEPEN + CMET*FU; // what the receptor sees
double SD = lep_signal(LEFF, LEC50, LHILL, SOCS, KSOCS, SREFB);
double DEFC = 1.0 - SD; // the deficit that drives hyperphagia
double OCC = LTOT/(LTOT + KDLEP); // receptor-occupancy readout
// ===========================================================================
// 2. INTAKE AND INFLUX
// ===========================================================================
double CGLP = 1000.0*GCEN/VGLP; // ng/mL
double EGLP = EMAXGLP*CGLP/(CGLP + EC50GLP);
double diet = (SOLVERTIME >= DIETT) ? 1.0 : 0.0;
double ffat = diet ? DIETFATN : FFAT;
double fcho = diet ? DIETCHON : FCHO;
// Two anorexigenic levers act on ONE intake term, so they are automatically
// less than additive (scenario 22).
double EITG = (EICLAMP > 0) ? EICLAMP : EI0*(1.0 + HMAX*DEFC)*(1.0 - EGLP);
double INSTOT = INS + INSX;
double insdr = INSTOT/(INSTOT + KINS50);
double IRSEL = fmax(0.0, IRX - 1.0)/(fmax(0.0, IRX - 1.0) + 3.0); // selective hepatic IR
double DFAT = EIS*ffat/9.0;
double CHOG = EIS*fcho/4.0;
double DNL = KDNL*CHOG*insdr*(1.0 + KDNLSEL*IRSEL)
*(1.0 + KLDNL*fmax(0.0, SREFB - SLAG));
// Adipose lipolytic release. This flux is RECYCLED: it is presented back to
// the buffer and must therefore appear on both sides of the balance, otherwise
// the futile cycle would look like free storage capacity.
double lipsup = 1.0/(1.0 + KINSN*INSTOT/8.0);
double JREL1 = KLIPSC*(1.0 + BASLIP)*ASC1*1000.0*lipsup;
double JREL2 = KLIPSC*(1.0 + BASLIP)*ASC2*1000.0*lipsup;
double JRELV = KLIPVI*(1.0 + BASLIP)*AVIS*1000.0*lipsup;
double JRELT = JREL1 + JREL2 + JRELV;
double JIN = DFAT + DNL + JRELT;
// ===========================================================================
// 3. THE BUFFER: CAPACITY-LIMITED ESTERIFICATION
// ===========================================================================
double ISENS = 1.0/fmax(IRX, IRMIN);
double jsat = JIN/(JIN + KMJ);
// Adipose esterification is driven by insulin ACTION = level x sensitivity,
// normalised to 1 in the healthy reference subject. This matters: in severe
// insulin resistance the adipocyte sees a 4-5x higher insulin concentration,
// so the signalling term is nearly preserved and CAPACITY -- not signalling --
// is what limits storage. Gating storage on sensitivity alone would empty the
// residual depots, which is the opposite of the hypertrophic adipocytes seen
// in partial lipodystrophy.
double xact = INSTOT*fmin(ISENS, 2.0);
double insact = 2.0*xact/(xact + KINSACT);
double hr1 = fmax(0.0, 1.0 - fill1) + HRFLOOR;
double hr2 = fmax(0.0, 1.0 - fill2) + HRFLOOR;
double hrv = fmax(0.0, 1.0 - fillv) + HRFLOOR;
double JST1 = KEST*C1*hr1*jsat*insact;
double JST2 = KEST*C2*hr2*jsat*insact;
double JSTV = KEST*CV*hrv*jsat*insact;
double JSTT = JST1 + JST2 + JSTV;
// Whole-body oxidation, with an adaptive component: patients with a large
// ectopic load up-regulate fat oxidation, which is why hyperphagia does not
// simply make them obese.
double BW = LEAN0 + FATM + (LTG + IMCL)/1000.0;
double OXW = KOXW*pow(BW, 0.75)*(1.0 + KOXAD*OXAD);
// THE overflow flux. Everything downstream of here is a consequence.
double JOV = fmax(0.0, JIN - JSTT - OXW);
// ===========================================================================
// 4. DRUG PK
// ===========================================================================
double CPIO = 1000.0*PCEN/VPIO; // ng/mL
double EPIO = CPIO/(CPIO + EC50PIO);
double CVOLT = VTIS/VVOLT; // mg/L
double CFIB = FCEN/VFIB; // mg/L
double EFIB = CFIB/(CFIB + EC50FIB);
double CEVI = EVC/VEVI; // mg/L
double CMETF = MTCEN/VMETF; // mg/L
double EMETF = EMAXMETF*CMETF/(CMETF + EC50METF);
double CSGL = SGCEN/VSGL; // mg/L
double UGE = UGEMAX*(CSGL/(CSGL + EC50SGL))*fmax(0.0, GLU - GSGLTHR)/GLU;
dxdt_MSC = -KAMET*MSC;
dxdt_MCEN = KAMET*MSC*FMET - (CLMET/VMETC)*MCEN - (QMET/VMETC)*MCEN + (QMET/VMETP)*MPER;
dxdt_MPER = (QMET/VMETC)*MCEN - (QMET/VMETP)*MPER;
dxdt_NAB = KNABON*NEUT*(CMET/(CMET + CNAB50))*(NABMAX - NAB) - KNABOFF*NAB;
dxdt_PGUT = -KAPIO*PGUT;
dxdt_PCEN = KAPIO*PGUT - (CLPIO/VPIO)*PCEN;
dxdt_VSC = -KAVOL*VSC;
dxdt_VCEN = KAVOL*VSC - (CLVOL/VVOLC)*VCEN - (QVOL/VVOLC)*VCEN + (QVOL/VVOLT)*VTIS;
dxdt_VTIS = (QVOL/VVOLC)*VCEN - (QVOL/VVOLT)*VTIS - KOUTVT*VTIS;
double kin3 = KOUT3*(1.0 - KFIB3*EFIB);
dxdt_APOC3 = kin3*(1.0 - EMAX3*CVOLT/(CVOLT + EC503)) - KOUT3*APOC3;
dxdt_FGUT = -KAFIB*FGUT;
dxdt_FCEN = KAFIB*FGUT - (CLFIB/VFIB)*FCEN;
dxdt_EVC = -(CLEVI/VEVI)*EVC;
dxdt_ANG3 = KOUTA*(1.0 - EMAXA*CEVI/(CEVI + EC50A)) - KOUTA*ANG3;
dxdt_OM3 = (((SOLVERTIME >= OM3T) ? OM3DOSE : 0.0) - OM3)/TAUOM3;
dxdt_MTGUT = -KAMET2*MTGUT;
dxdt_MTCEN = KAMET2*MTGUT - (CLMETF/VMETF)*MTCEN;
double insinf = (SOLVERTIME >= INST) ? INSDOSE : 0.0;
dxdt_ISCD = insinf - KAI*ISCD;
dxdt_INSX = KAI*ISCD*SCALEI - KINSCL*INSX; // exogenous: no hepatic first pass
dxdt_GSC = -KAGLP*GSC;
dxdt_GCEN = KAGLP*GSC - (CLGLP/VGLP)*GCEN;
dxdt_SGGUT = -KASGL*SGGUT;
dxdt_SGCEN = KASGL*SGGUT - (CLSGL/VSGL)*SGCEN;
// ===========================================================================
// 5. LEPTIN-AXIS STATES
// ===========================================================================
dxdt_LEPEN = (lep_ss - LEPEN)/TAULEP;
dxdt_SOCS = (SD - SOCS)/TAUSOCS; // the brake tracks the sustained signal
dxdt_SLAG = (SD - SLAG)/TAUSLAG;
dxdt_EIS = (EITG - EIS)/TAUEI;
// ===========================================================================
// 6. ADIPOSE DEPOTS, CAPACITY, ADAPTIVE OXIDATION
// ===========================================================================
dxdt_ASC1 = (JST1 - JREL1)/1000.0;
dxdt_ASC2 = (JST2 - JREL2)/1000.0;
dxdt_AVIS = (JSTV - JRELV)/1000.0;
dxdt_CAPM = ((1.0 + PREAD*KCAPM*EPIO) - CAPM)/TAUCAP;
dxdt_OXAD = ((JOV/JOVREF)/(1.0 + JOV/JOVREF) - OXAD)/TAUOXAD;
// ===========================================================================
// 7. ECTOPIC LIPID (the overflow cascade)
// ===========================================================================
// VLDL-TG secretion is written as PURELY SUBSTRATE-DRIVEN (Adiels): there is no
// direct leptin term on it. An earlier draft gave leptin a primary inhibitory
// action here, and that draft predicted that leptin replacement at CLAMPED food
// intake makes steatosis 71% WORSE -- because suppressing the largest hepatic
// disposal route at unchanged substrate delivery has to raise the pool. Since
// the observed effect of leptin at fixed intake is the opposite, the falling
// VLDL secretion seen on metreleptin must be a CONSEQUENCE of the smaller
// hepatic triglyceride pool, not a cause of it. KLVLDL is retained at 0 so the
// discarded structure stays visible and testable.
double vldl = VMVLDL*(LTG/(LTG + KMLTG))
*(1.0 + KLVLDL*fmax(0.0, SREFB - SLAG))
*(1.0 - KOM3VLDL*OM3)
*(1.0 - KEVVLDL*(1.0 - ANG3));
// Leptin acts on the liver in three directions at once (AMPK-driven FA
// oxidation up, SREBP-1c-driven lipogenesis down, VLDL secretion down). Only
// modelling the third would predict that leptin replacement at fixed food
// intake makes steatosis WORSE, which is not what pair-fed animals or the
// Petersen hepatic-TG data show.
double hepox = VMOX*(1.0 + KFIBOX*EFIB + KOM3OX*OM3)
*(1.0 - KLOX*fmax(0.0, SREFB - SLAG))*LTG/(LTG + KMOX);
// The hepatic SHARE of overflow is not a constant: as adipose disappears the
// liver becomes a larger fraction of the remaining sink.
// The liver competes with SUBCUTANEOUS adipose only. Visceral fat is not a
// competing sink at all -- it drains into the portal vein, so it delivers TO
// the liver. This is why partial lipodystrophy, which trades subcutaneous for
// visceral fat, is steatotic despite a near-normal total fat mass.
double scfr0 = fmin(1.0, (ASC1 + ASC2)/SCFATN);
double flive = fmin(0.88, FLIV + KFLIVA*(1.0 - scfr0));
// Visceral lipolysis drains into the portal vein, so it reaches the liver at a
// higher fraction than systemic NEFA does. This is the term that makes
// PARTIAL lipodystrophy (visceral expansion, subcutaneous loss) steatotic.
dxdt_LTG = FHNEFA*(JREL1 + JREL2) + FHVIS*JRELV + flive*JOV + DNL*0.35
- hepox - vldl;
dxdt_IMCL = FMUS*JOV - KIMOX*IMCL;
dxdt_PLIP = ((FPANC*JOV)/(FPANC*JOVREF) - PLIP)/TAUECT;
dxdt_MYOL = ((FMYO*JOV)/(FMYO*JOVREF) - MYOL)/TAUECT;
dxdt_KLIP = ((FKID*JOV)/(FKID*JOVREF) - KLIP)/TAUECT;
// ===========================================================================
// 8. PLASMA LIPID
// ===========================================================================
double fatfrac = fmin(1.5, FATM/FATN);
double lplact = (LPLFLOOR + (1.0 - LPLFLOOR)*fatfrac)
*(1.0 + KAP3*(1.0 - APOC3))
*(1.0 + KANG*(1.0 - ANG3))
*(1.0 + KFIBL*EFIB);
double tgin = (vldl + FPL*JOV)*1000.0/VPDL;
double tgout = VMTG*lplact*PTG/(KMTGP + PTG) + KLINTG*PTG;
dxdt_PTG = tgin - tgout;
double nefa_ss = NEFA0*(1.0 + KLIPB*BASLIP)*lipsup*(1.0 + 0.5*fmax(0.0, JOV/JOVREF - 1.0));
dxdt_NEFA = (nefa_ss - NEFA)/TAUNEFA;
double adpn_ss = ADPNFLR + ADPN0*fatfrac/(1.0 + 0.5*fmax(0.0, IRX - 1.0))
+ KADPPIO*EPIO*fmin(1.0, fatfrac + 0.15);
dxdt_ADPN = (adpn_ss - ADPN)/TAUADP;
// ===========================================================================
// 9. INSULIN RESISTANCE, GLUCOSE, BETA CELL
// ===========================================================================
double IMCLNN = IMCL/IMCLN0;
double IRnow = 1.0 + AHFF*(HFF - 3.0) + BIM*(IMCLNN - 1.0) + CINFL*(INFL - 0.17)
- CLEP*(SLAG - SDREF) - DADP*(ADPN - ADPN0)/ADPN0;
IRnow = fmax(IRMIN, IRnow);
dxdt_IRX = (IRnow - IRX)/TAUIR;
// SELECTIVE hepatic insulin resistance, written with the SAME term that keeps
// lipogenesis insulin-driven (IRSEL). One coefficient therefore produces both
// halves of the paradox: glucose production escapes insulin while de novo
// lipogenesis still obeys it. This is why the model can be hyperinsulinaemic
// AND hyperglycaemic AND hyperlipogenic at once.
double hgp = HGP0*(1.0 + KHIR*fmax(0.0, IRX - 1.0))
/(1.0 + KHINS*INSTOT*(1.0 - IRSEL))*(1.0 - EMETF);
double rd = (KGZ + KGI*INSTOT*fmin(ISENS, 2.0))*GLU;
double uger = KREN*fmax(0.0, GLU - GRENTHR);
dxdt_GLU = hgp - rd - uger - UGE;
double isec = BCF*SMAX*(pow(GLU,2.0)/(pow(G50,2.0) + pow(GLU,2.0)))
*(1.0 + KCOMP*fmax(0.0, IRX - 1.0));
// Hepatic insulin extraction falls as insulin resistance rises, which is part
// of why lipodystrophy is hyperinsulinaemic out of proportion to secretion.
double kinscl = KINSCL/(1.0 + KIEXT*fmax(0.0, IRX - 1.0));
dxdt_INS = isec - kinscl*INS;
dxdt_BCF = KBREC*(1.0 - BCF)
- BCF*(KBLIP*fmax(0.0, PLIP - 1.0)
+ KBGLU*fmax(0.0, (GLU - GTOXTHR)/100.0));
dxdt_A1C = (((MGFAC*GLU) + 46.7)/28.7 - A1C)/TAUA1C;
// ===========================================================================
// 10. LIVER INJURY, FIBROSIS, KIDNEY, ANDROGENS, PANCREATITIS
// ===========================================================================
double alt_ss = ALT0 + KALT*fmax(0.0, HFF - 3.0)*(1.0 + 0.5*INFL);
dxdt_ALT = (alt_ss - ALT)/TAUALT;
dxdt_INFL = (HFF/(HFF + K50INFL) - INFL)/TAUINFL;
dxdt_FIB = KFIB*fmax(0.0, INFL - INFLTHR)*(1.0 - FIB/4.0) - KFREG*FIB;
double prot_ss = 0.05 + KPR*fmax(0.0, IRX - 1.0) + KPRG*fmax(0.0, (GLU - 140.0)/100.0)
+ KPRL*fmax(0.0, KLIP - 1.0);
dxdt_PROT = (prot_ss - PROT)/TAUPROT;
dxdt_EGFR = -KEGFR*PROT*(EGFR > 10.0 ? 1.0 : 0.0);
double and_ss = 1.0 + KAND*fmax(0.0, (INSTOT - 8.0)/8.0)/(1.0 + fmax(0.0, (INSTOT - 8.0)/8.0)/4.0);
dxdt_ANDX = (and_ss - ANDX)/TAUAND;
dxdt_TAT = 1.0/(1.0 + exp(-(PTG - TGTHR)/50.0));
double tgr = pow(PTG/1000.0, 4.0);
dxdt_PANCH = KPANC*tgr/(1.0 + tgr);
$TABLE
// Recomputed from the states at the output time so that captured values never
// depend on where the solver happened to stop.
double t_fill = fmin(1.0, (ASC1 + ASC2)/(CAP1N*FC1*CAPM + CAP2N*FC2*CAPM + 1e-9));
double t_CMET = MCEN/VMETC;
double t_FU = 1.0/(1.0 + NAB);
double t_LTOT = LEPEN + t_CMET;
double t_LEFF = LEPEN + t_CMET*t_FU;
double t_SD = lep_signal(t_LEFF, LEC50, LHILL, SOCS, KSOCS, SREFB);
double t_DEF = 1.0 - t_SD;
double t_OCC = t_LTOT/(t_LTOT + KDLEP);
double t_HFF = HFF;
double t_FAT = FATM;
double t_INST = INS + INSX;
double t_ISEN = 1.0/fmax(IRX, IRMIN);
double t_diet = (TIME >= DIETT) ? 1.0 : 0.0;
double t_ffat = t_diet ? DIETFAT : FFAT;
double t_fcho = t_diet ? (FFAT + FCHO - DIETFAT) : FCHO;
double t_insdr = t_INST/(t_INST + KINS50);
double t_IRSEL = fmax(0.0, IRX - 1.0)/(fmax(0.0, IRX - 1.0) + 3.0);
double t_DNL = KDNL*(EIS*t_fcho/4.0)*t_insdr*(1.0 + KDNLSEL*t_IRSEL)
*(1.0 + KLDNL*fmax(0.0, SREFB - SLAG));
double t_lsup = 1.0/(1.0 + KINSN*t_INST/8.0);
double t_JREL = (KLIPSC*(ASC1 + ASC2) + KLIPVI*AVIS)*(1.0 + BASLIP)*1000.0*t_lsup;
double t_JIN = EIS*t_ffat/9.0 + t_DNL + t_JREL;
double t_jsat = t_JIN/(t_JIN + KMJ);
double t_xact = t_INST*fmin(t_ISEN, 2.0);
double t_iact = 2.0*t_xact/(t_xact + KINSACT);
double t_C1 = CAP1N*FC1*CAPM;
double t_C2 = CAP2N*FC2*CAPM;
double t_CV = CAPVN*FCV;
double t_JST = KEST*t_jsat*t_iact*
( t_C1*(fmax(0.0, 1.0 - ASC1/(t_C1 + 1e-9)) + HRFLOOR)
+ t_C2*(fmax(0.0, 1.0 - ASC2/(t_C2 + 1e-9)) + HRFLOOR)
+ t_CV*(fmax(0.0, 1.0 - AVIS/(t_CV + 1e-9)) + HRFLOOR) );
double t_BW = LEAN0 + FATM + (LTG + IMCL)/1000.0;
double t_OXW = KOXW*pow(t_BW, 0.75)*(1.0 + KOXAD*OXAD);
double t_JOV = fmax(0.0, t_JIN - t_JST - t_OXW);
double t_CAP = t_C1 + t_C2 + t_CV;
double t_fatf = fmin(1.0, (ASC1 + ASC2)/SCFATN);
double t_fliv = fmin(0.88, FLIV + KFLIVA*(1.0 - t_fatf));
double t_JHEP = FHNEFA*(KLIPSC*(ASC1 + ASC2)*(1.0 + BASLIP)*1000.0*t_lsup)
+ FHVIS*(KLIPVI*AVIS*(1.0 + BASLIP)*1000.0*t_lsup)
+ t_fliv*t_JOV + 0.35*t_DNL;
$CAPTURE @annotated
t_HFF : hepatic fat fraction (%)
t_FAT : total body fat (kg)
t_CAP : total storage capacity (kg)
t_LTOT : total leptin, assay-equivalent (ng/mL)
t_LEFF : receptor-available leptin (ng/mL)
t_SD : hypothalamic leptin signal (0-1)
t_DEF : signal deficit = 1 - SD
t_OCC : leptin receptor occupancy readout (0-1)
t_JIN : lipid influx presented to the buffer (g/d)
t_JST : total adipose storage flux (g/d)
t_JOV : OVERFLOW flux (g/d)
t_OXW : whole-body fat oxidation (g/d)
t_BW : body weight (kg)
t_INST : total insulin, endogenous + exogenous (uU/mL)
t_ISEN : insulin sensitivity (1/IR)
t_fill : mean subcutaneous fractional fill
t_DNL : de novo lipogenesis (g/d)
t_JREL : adipose lipolytic release, recycled (g/d)
t_JHEP : fatty-acid flux delivered to the liver (g/d)
t_fliv : hepatic share of the overflow
$OMEGA @annotated @block
ELEAN : 0 : between-subject variability on maintenance intake
$SIGMA 0
)---"
## -----------------------------------------------------------------------------
## BUILD
## -----------------------------------------------------------------------------
mod <- mcode_cache("lipo_qsp", lipo_code, atol = 1e-8, rtol = 1e-8, maxsteps = 500000)
## =============================================================================
## PHENOTYPES
##
## A phenotype in this model is NOT a label. It is a capacity distribution
## (FC1/FC2/FCV), an expandability reserve (PREAD) and a secretion coefficient
## (LEPKM). Everything else -- leptin level, hepatic fat, triglyceride,
## HbA1c -- is an OUTPUT of running those numbers to steady state.
##
## The two control phenotypes matter as much as the disease ones:
## lep_def fat mass HIGH, signal ZERO -> the pure signal-deficit arm
## obese fat mass HIGH, signal HIGH -> the pure capacity-exceeded arm
## Together they show that the model's two deficits are separable.
## =============================================================================
phenotypes <- list(
normal = list(FC1 = 1.00, FC2 = 1.00, FCV = 1.00, PREAD = 1.00, LEPKM = 1.0,
LEAN0 = 45, label = "Healthy control"),
cgl = list(FC1 = 0.02, FC2 = 0.05, FCV = 0.05, PREAD = 0.03, LEPKM = 1.0,
LEAN0 = 48, label = "CGL (AGPAT2/BSCL2)"),
agl = list(FC1 = 0.05, FC2 = 0.08, FCV = 0.15, PREAD = 0.15, LEPKM = 1.0,
LEAN0 = 46, label = "Acquired generalised (Lawrence)"),
fpld_low = list(FC1 = 0.12, FC2 = 0.90, FCV = 1.70, PREAD = 0.55, LEPKM = 1.0,
LEAN0 = 45, label = "FPLD2, low leptin"),
fpld_high = list(FC1 = 0.22, FC2 = 1.30, FCV = 1.80, PREAD = 0.60, LEPKM = 1.0,
LEAN0 = 45, label = "FPLD2, near-normal leptin"),
apl = list(FC1 = 0.90, FC2 = 0.15, FCV = 0.60, PREAD = 0.50, LEPKM = 1.0,
LEAN0 = 45, label = "Acquired partial (Barraquer-Simons)"),
plin1 = list(FC1 = 0.30, FC2 = 0.55, FCV = 1.30, PREAD = 0.45, LEPKM = 1.0,
BASLIP = 0.8, LEAN0 = 45, label = "PLIN1 FPLD4 (high basal lipolysis)"),
lep_def = list(FC1 = 1.60, FC2 = 1.60, FCV = 1.60, PREAD = 1.00, LEPKM = 0.02,
LEAN0 = 45, label = "Congenital leptin deficiency (control)"),
obese = list(FC1 = 1.80, FC2 = 1.80, FCV = 2.00, PREAD = 1.00, LEPKM = 1.0,
LEAN0 = 55, label = "Common obesity (control)")
)
pheno_param <- function(p) {
p$label <- NULL
p
}
## =============================================================================
## BASELINE: each phenotype is run to ITS OWN steady state
##
## Hand-typed initial conditions would let a phenotype be asserted rather than
## derived, and would make every treatment contrast depend on how good the
## guess was. Instead every scenario starts from a 12-year burn-in of the
## untreated phenotype.
## =============================================================================
## Run an arbitrary parameter set to its own steady state. The Shiny app uses
## this so that moving a capacity slider produces a genuinely re-derived
## baseline rather than a disease phenotype with one number edited.
baseline_pheno <- function(p, burn = 12*365) {
m <- param(mod, p)
out <- mrgsim_df(zero_re(m), end = burn, delta = burn, recsort = 3)
last <- out[nrow(out), ]
st <- as.list(last[intersect(names(last), names(as.list(init(mod))))])
list(mod = init(m, st), init = st, param = p)
}
baseline <- function(pheno = "cgl", burn = 12*365) {
p <- pheno_param(phenotypes[[pheno]])
m <- param(mod, p)
out <- mrgsim_df(zero_re(m), end = burn, delta = burn, recsort = 3)
last <- out[nrow(out), ]
st <- as.list(last[intersect(names(last), names(as.list(init(mod))))])
list(mod = init(m, st), init = st, param = p)
}
baseline_table <- function() {
rows <- lapply(names(phenotypes), function(nm) {
b <- baseline(nm)
s <- mrgsim_df(zero_re(b$mod), end = 1, delta = 1)
s <- s[nrow(s), ]
data.frame(phenotype = phenotypes[[nm]]$label,
fat_kg = round(s$t_FAT, 2),
capacity_kg = round(s$t_CAP, 2),
leptin = round(s$t_LTOT, 2),
signal_SD = round(s$t_SD, 3),
deficit = round(s$t_DEF, 3),
intake_kcal = round(s$EIS),
JIN_g = round(s$t_JIN, 1),
JOV_g = round(s$t_JOV, 1),
HFF_pct = round(s$t_HFF, 1),
TG = round(s$PTG),
glucose = round(s$GLU),
HbA1c = round(s$A1C, 2),
insulin = round(s$INS, 1),
IR = round(s$IRX, 2),
betacell = round(s$BCF, 2),
adipo = round(s$ADPN, 1),
ALT = round(s$ALT),
fibrosis = round(s$FIB, 2),
protein = round(s$PROT, 2),
row.names = NULL)
})
do.call(rbind, rows)
}
## =============================================================================
## DOSING HELPERS
##
## Metreleptin amounts are in micrograms so that amount/volume comes out in
## ng/mL, the unit every leptin paper reports.
## =============================================================================
met_sc <- function(mgkg = 0.06, wt = 60, start = 730, dur = 3650)
ev(amt = mgkg*wt*1000, cmt = "MSC", time = start, ii = 1, addl = dur - 1)
met_fixed <- function(mg = 5, start = 730, dur = 3650)
ev(amt = mg*1000, cmt = "MSC", time = start, ii = 1, addl = dur - 1)
pio_oral <- function(mg = 45, start = 730, dur = 3650)
ev(amt = mg, cmt = "PGUT", time = start, ii = 1, addl = dur - 1)
vol_sc <- function(mg = 300, start = 730, dur = 3650)
ev(amt = mg, cmt = "VSC", time = start, ii = 7, addl = floor(dur/7) - 1)
fen_oral <- function(mg = 145, start = 730, dur = 3650)
ev(amt = mg, cmt = "FGUT", time = start, ii = 1, addl = dur - 1)
evi_iv <- function(mg = 900, start = 730, n = 130)
ev(amt = mg, cmt = "EVC", time = start, ii = 28, addl = n - 1)
metf_oral <- function(mg = 1000, start = 730, dur = 3650)
ev(amt = mg, cmt = "MTGUT", time = start, ii = 0.5, addl = dur*2 - 1)
glp_sc <- function(mg = 1, start = 730, dur = 3650)
ev(amt = mg, cmt = "GSC", time = start, ii = 7, addl = floor(dur/7) - 1)
sgl_oral <- function(mg = 10, start = 730, dur = 3650)
ev(amt = mg, cmt = "SGGUT", time = start, ii = 1, addl = dur - 1)
## =============================================================================
## SCENARIO ENGINE
##
## Every scenario is the same 53 equations started from the same phenotype
## steady state. Contrasts are therefore differences in the intervention
## only -- which is what makes 7 vs 8, 13 vs 14 and 6 vs 11 interpretable.
## =============================================================================
run_scn <- function(label, pheno, events = NULL, pset = NULL,
end = 12*365, delta = 7) {
b <- baseline(pheno)
m <- zero_re(b$mod)
if (!is.null(pset)) m <- param(m, pset)
out <- if (is.null(events)) mrgsim_df(m, end = end, delta = delta)
else mrgsim_df(m, events = events, end = end, delta = delta)
out$scenario <- label
out$pheno <- pheno
out
}
scenarios <- list(
## ---- 1-5: natural history and the two control phenotypes ---------------
s1 = function() run_scn("1 CGL untreated", "cgl"),
s2 = function() run_scn("2 FPLD2 low-leptin untreated", "fpld_low"),
s3 = function() run_scn("3 Healthy control", "normal"),
s4 = function() run_scn("4 Common obesity (capacity exceeded, signal high)", "obese"),
s5 = function() run_scn("5 Congenital leptin deficiency (signal zero, capacity intact)", "lep_def"),
## ---- 6-9: the metreleptin dose/deficit dissociation --------------------
s6 = function() run_scn("6 CGL + metreleptin 0.06 mg/kg", "cgl",
met_sc(0.06, 60, start = 730)),
s7 = function() run_scn("7 FPLD low-leptin + metreleptin 0.06", "fpld_low",
met_sc(0.06, 60, start = 730)),
s8 = function() run_scn("8 FPLD near-normal leptin + metreleptin 0.06", "fpld_high",
met_sc(0.06, 60, start = 730)),
s9 = function() run_scn("9 FPLD near-normal leptin + metreleptin 0.13 (dose escalation)",
"fpld_high", met_sc(0.13, 60, start = 730)),
## ---- 10-12: reversibility, immunogenicity, pair-feeding ---------------
s10 = function() run_scn("10 CGL + metreleptin, withdrawn at year 6", "cgl",
met_sc(0.06, 60, start = 730, dur = 4*365)),
s11 = function() run_scn("11 CGL + metreleptin with neutralising ADA", "cgl",
met_sc(0.06, 60, start = 730), pset = list(NEUT = 1)),
s12 = function() {
## Pair-fed clamp: intake is held at the untreated steady-state value, so
## whatever effect remains is the DIRECT peripheral action of leptin.
b <- baseline("cgl")
ei <- b$init$EIS
run_scn("12 CGL + metreleptin, food intake clamped (pair-fed)", "cgl",
met_sc(0.06, 60, start = 730), pset = list(EICLAMP = ei))
},
s12b = function() {
b <- baseline("cgl")
ei <- b$init$EIS
run_scn("12b CGL untreated, intake clamped (pair-fed control)", "cgl",
NULL, pset = list(EICLAMP = ei))
},
## ---- 13-14: capacity-directed therapy needs a substrate --------------
s13 = function() run_scn("13 CGL + pioglitazone 45 mg (no preadipocyte substrate)",
"cgl", pio_oral(45, start = 730)),
s14 = function() run_scn("14 FPLD low-leptin + pioglitazone 45 mg (substrate present)",
"fpld_low", pio_oral(45, start = 730)),
## ---- 15-17: the spillover agents ------------------------------------