qsp/fabry-disease/fbr_mrgsolve_model.R at main · pipetcpt/qsp

Repository navigation

History

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

##############################################################################

# Fabry Disease QSP Model — mrgsolve ODE (R)

# Disease : Fabry Disease (파브리병, α-Galactosidase A Deficiency)

# Drug targets: ERT (agalsidase alfa/beta, pegunigalsidase alfa),

# Migalastat (pharmacological chaperone, amenable mutations),

# Lucerastat (GCS inhibitor, substrate reduction therapy)

# Compartments: 22 ODE compartments

# Scenarios: 6 treatment scenarios + natural history

# Parameter calibration: Fabry Registry, FABRY-001 trial, ATTRACT trial,

# BRIGHT trial, MODIFY trial (lucerastat)

# Author: Claude Code Routine (CCR) | Date: 2026-06-24

##############################################################################

library(mrgsolve)

library(dplyr)

library(ggplot2)

library(tidyr)

##############################################################################

# MODEL DEFINITION

##############################################################################

fbr_code <- '

$PROB

Fabry Disease QSP Model

- X-linked alpha-galactosidase A (GLA) deficiency

- Gb3/lyso-Gb3 multi-organ accumulation

- Drug treatments: ERT (agalsidase beta/alfa/pegunigalsidase alfa),

Migalastat (chaperone, amenable mutations),

Lucerastat (GCS inhibitor, SRT)

- 22 ODE compartments

- Reference: Schiffmann R 2001 Ann Intern Med; Mehta A 2009 Eur J Clin Invest;

Hughes DA 2017 Lancet; Germain DP 2016 NEJM (ATTRACT);

Linthorst GE 2010 J Inherit Metab Dis

$PARAM

// ── ERT (Agalsidase Beta: Fabrazyme 1 mg/kg Q2W IV) ──────────────────────

CL_AGAB = 0.42 // L/h clearance agalsidase beta (BW-scaled)

Vc_AGAB = 2.8 // L central volume

Vp_AGAB = 8.1 // L peripheral volume

Q_AGAB = 1.25 // L/h intercompartmental clearance

k_lys_AGAB = 0.35 // 1/h M6PR-mediated lysosomal uptake from plasma

CL_lys_AGAB = 0.18 // 1/h lysosomal ERT degradation

// Reference: Vedder 2007 J Am Soc Nephrol; t1/2 plasma ~45 min

// ── ERT (Agalsidase Alfa: Replagal 0.2 mg/kg Q2W IV) ────────────────────

CL_AGAA = 0.55 // L/h (higher CL than agaB, lower dose)

Vc_AGAA = 3.1 // L

Vp_AGAA = 7.8 // L

Q_AGAA = 1.10 // L/h

k_lys_AGAA = 0.28 // 1/h

CL_lys_AGAA = 0.20 // 1/h

// Reference: Schiffmann 2001 Ann Intern Med

// ── Migalastat PK (150 mg oral QOD) ──────────────────────────────────────

ka_MIG = 0.82 // 1/h absorption rate constant

F_MIG = 0.75 // bioavailability

Vc_MIG = 48.0 // L/kg central volume (t1/2 ~3.5h)

CL_MIG = 9.5 // L/h total clearance

// Reference: Germain 2016 NEJM; Bichet 2017 J Inherit Metab Dis

// ── Lucerastat PK (SRT, 1000 mg TID oral) ────────────────────────────────

ka_LUC = 0.60 // 1/h absorption

F_LUC = 0.65 // bioavailability

CL_LUC = 12.5 // L/h

Vc_LUC = 95.0 // L (t1/2 ~8h)

IC50_LUC = 0.18 // μg/mL GCS inhibition IC50

Emax_LUC = 0.42 // maximum fractional GCS inhibition

// ── α-Gal A Enzyme Biology ───────────────────────────────────────────────

E_GalA_base = 8.0 // nmol/h/mg normal enzyme activity (hemizygote male = 0)

E_GalA0_classic = 0.05 // nmol/h/mg classic Fabry (residual <1%)

E_GalA0_late = 2.5 // nmol/h/mg late-onset Fabry (residual ~30%)

Km_ERT = 1.5 // ng/mL Michaelis constant for ERT→enzyme activity

Emax_ERT = 72.0 // nmol/h/mg max enzyme activity restored by ERT

// Migalastat PD

Emax_MIG = 6.0 // nmol/h/mg max enzyme activity restored (amenable mut.)

EC50_MIG = 0.25 // μg/mL EC50 migalastat

// ── Gb3 / Lyso-Gb3 Dynamics ──────────────────────────────────────────────

// Plasma Gb3

kin_Gb3 = 0.35 // μg/mL/h Gb3 synthesis rate (healthy ~0 accumulation)

kel_Gb3_norm = 0.18 // 1/h normal Gb3 elimination (enzyme-dependent)

kel_Gb3_base = 0.008 // 1/h residual Gb3 elimination (passive)

Gb3_PL_ss = 0.8 // μg/mL Fabry patient plasma Gb3 steady-state (~healthy <0.1)

// Lyso-Gb3 plasma

kin_LGB3 = 0.025 // μg/mL/h lyso-Gb3 generation from Gb3

kel_LGB3 = 0.15 // 1/h lyso-Gb3 clearance

LGB3_ss = 12.0 // μg/L classic male Fabry patient lyso-Gb3 at SS (~healthy <2)

// ── Kidney Gb3 / eGFR ────────────────────────────────────────────────────

kin_Gb3k = 0.010 // nmol/mg Cr/h kidney Gb3 accumulation rate

kel_Gb3k = 0.004 // 1/h kidney Gb3 clearance (enzyme-dependent)

Gb3k_ss0 = 50.0 // nmol/mg Cr classic Fabry kidney Gb3 at baseline

eGFR0 = 90.0 // mL/min/1.73m² baseline eGFR (young adult)

k_eGFR_Gb3 = 0.0003 // fractional eGFR decline per unit Gb3k per year

UPCR0 = 350.0 // mg/g Cr baseline UPCR in advanced Fabry nephropathy

// ── Cardiac Gb3 / LVMi ───────────────────────────────────────────────────

kin_Gb3h = 0.008 // μg/mg Cr equivalent/h cardiac Gb3 accumulation

kel_Gb3h = 0.003 // 1/h

Gb3h_ss0 = 45.0 // μg/mg (normalized units) baseline cardiac Gb3

LVMi0 = 148.0 // g/m² baseline LVMi in symptomatic Fabry cardiomyopathy

k_LVMi_Gb3h = 0.0008 // LVMi change per unit cardiac Gb3

LVMi_min = 95.0 // g/m² minimum achievable (treatment target)

// ── Neurological / Neuropathic Pain ──────────────────────────────────────

Pain0 = 6.5 // BPI-SF score baseline (0–10)

k_pain_ERT = 0.25 // Pain reduction rate constant with effective ERT

k_pain_base = 0.003 // spontaneous pain increase rate

// ── Inflammatory Index ────────────────────────────────────────────────────

Inflam0 = 1.8 // arbitrary units elevated in Fabry (IL-6/TNF-alpha)

k_inflam_in = 0.020 // inflammatory driving from Gb3

k_inflam_out = 0.18 // inflammatory resolution rate

// ── GCS Inhibition Chain ──────────────────────────────────────────────────

GCS_base_act = 1.0 // normalized GCS activity

k_Gb3_GCS = 0.95 // fraction of Gb3 synthesis through GCS pathway

// ── Disease phenotype ────────────────────────────────────────────────────

IS_CLASSIC = 1 // 1=classic Fabry, 0=late-onset

IS_AMENABLE = 0 // 1=amenable mutation (migalastat), 0=non-amenable

BW = 70 // kg body weight

$CMT

// ERT Compartments (Agalsidase Beta, dominant ERT in model)

A_AGAB_C // agalsidase beta: central plasma (μg)

A_AGAB_P // agalsidase beta: peripheral tissue (μg)

A_AGAB_LYS // agalsidase beta: lysosomal (active enzyme pool, μg)

// Agalsidase Alfa Compartments

A_AGAA_C // agalsidase alfa: central (μg)

A_AGAA_P // agalsidase alfa: peripheral (μg)

A_AGAA_LYS // agalsidase alfa: lysosomal (μg)

// Migalastat Compartments

A_MIG_GUT // migalastat: gut (μg)

A_MIG_C // migalastat: central plasma (μg)

// Lucerastat/SRT Compartments

A_LUC_GUT // lucerastat: gut (μg)

A_LUC_C // lucerastat: central plasma (μg)

// Enzyme Activity

E_GalA // α-galactosidase A activity (nmol/h/mg)

// Glycosphingolipids

GB3_PLM // plasma Gb3 (μg/mL)

GB3_KID // kidney Gb3 (nmol/mg Cr, normalized)

GB3_HRT // cardiac Gb3 (normalized units)

LGB3_PLM // plasma lyso-Gb3 (μg/L)

// Inflammatory Index

INFLAM // inflammation (normalized, IL-6/TNF-α composite)

// Organ Function

eGFR // glomerular filtration rate (mL/min/1.73m²)

UPCR // urine protein:creatinine ratio (mg/g)

LVMi // left ventricular mass index (g/m²)

PAIN // BPI-SF neuropathic pain score (0–10)

$MAIN

// ── Set baseline enzyme activity based on phenotype ─────────────────────

double E_GalA_init = IS_CLASSIC > 0.5 ? E_GalA0_classic : E_GalA0_late;

// Initial concentrations

E_GalA_0 = E_GalA_init;

GB3_PLM_0 = Gb3_PL_ss;

GB3_KID_0 = Gb3k_ss0;

GB3_HRT_0 = Gb3h_ss0;

LGB3_PLM_0 = LGB3_ss;

INFLAM_0 = Inflam0;

eGFR_0 = eGFR0;

UPCR_0 = UPCR0;

LVMi_0 = LVMi0;

PAIN_0 = Pain0;

$ODE

// ── ERT (Agalsidase Beta) PK ─────────────────────────────────────────────

double C_AGAB = A_AGAB_C / (Vc_AGAB * BW); // ng/mL

dxdt_A_AGAB_C = -CL_AGAB * C_AGAB * BW

- Q_AGAB * (C_AGAB - A_AGAB_P/(Vp_AGAB * BW))

- k_lys_AGAB * A_AGAB_C;

dxdt_A_AGAB_P = Q_AGAB * (C_AGAB * BW - A_AGAB_P/Vp_AGAB)

- A_AGAB_P * 0.05; // slow tissue clearance

dxdt_A_AGAB_LYS = k_lys_AGAB * A_AGAB_C

- CL_lys_AGAB * A_AGAB_LYS;

// ── ERT (Agalsidase Alfa) PK ──────────────────────────────────────────────

double C_AGAA = A_AGAA_C / (Vc_AGAA * BW);

dxdt_A_AGAA_C = -CL_AGAA * C_AGAA * BW

- Q_AGAA * (C_AGAA - A_AGAA_P/(Vp_AGAA * BW))

- k_lys_AGAA * A_AGAA_C;

dxdt_A_AGAA_P = Q_AGAA * (C_AGAA * BW - A_AGAA_P/Vp_AGAA)

- A_AGAA_P * 0.05;

dxdt_A_AGAA_LYS = k_lys_AGAA * A_AGAA_C - CL_lys_AGAA * A_AGAA_LYS;

// ── Migalastat PK ─────────────────────────────────────────────────────────

dxdt_A_MIG_GUT = -ka_MIG * A_MIG_GUT;

dxdt_A_MIG_C = ka_MIG * F_MIG * A_MIG_GUT

- (CL_MIG / Vc_MIG) * A_MIG_C;

double C_MIG = A_MIG_C / Vc_MIG; // μg/mL

// ── Lucerastat PK ─────────────────────────────────────────────────────────

dxdt_A_LUC_GUT = -ka_LUC * A_LUC_GUT;

dxdt_A_LUC_C = ka_LUC * F_LUC * A_LUC_GUT

- (CL_LUC / Vc_LUC) * A_LUC_C;

double C_LUC = A_LUC_C / Vc_LUC; // μg/mL

// ── α-Gal A Enzyme Activity ───────────────────────────────────────────────

// Combined ERT contribution from both agalsidase beta and alfa

double LYS_total = A_AGAB_LYS + A_AGAA_LYS; // μg in lysosomal pool

double E_ERT = Emax_ERT * LYS_total / (Km_ERT * BW + LYS_total);

// Migalastat contribution (only for amenable mutations)

double E_MIG = IS_AMENABLE > 0.5 ?

Emax_MIG * C_MIG / (EC50_MIG + C_MIG) : 0.0;

double E_GalA_target = E_GalA_init + E_ERT + E_MIG;

// Enzyme activity dynamics (turnover model)

dxdt_E_GalA = 0.15 * (E_GalA_target - E_GalA);

// ── GCS Inhibition (Lucerastat SRT) ──────────────────────────────────────

double GCS_inhib = Emax_LUC * C_LUC / (IC50_LUC + C_LUC);

double Gb3_syn_rate = kin_Gb3 * (1.0 - k_Gb3_GCS * GCS_inhib);

// ── Plasma Gb3 ───────────────────────────────────────────────────────────

double Gb3_elim_enz = kel_Gb3_norm * (E_GalA / (E_GalA + 1.0));

double Gb3_elim_total = kel_Gb3_base + Gb3_elim_enz;

dxdt_GB3_PLM = Gb3_syn_rate - Gb3_elim_total * GB3_PLM;

// ── Kidney Gb3 ────────────────────────────────────────────────────────────

double Gb3k_elim_enz = kel_Gb3k * (E_GalA / (E_GalA + 0.5));

dxdt_GB3_KID = kin_Gb3k - (0.002 + Gb3k_elim_enz) * GB3_KID;

// ── Cardiac Gb3 ──────────────────────────────────────────────────────────

double Gb3h_elim_enz = kel_Gb3h * (E_GalA / (E_GalA + 0.5));

dxdt_GB3_HRT = kin_Gb3h - (0.001 + Gb3h_elim_enz) * GB3_HRT;

// ── Plasma Lyso-Gb3 ──────────────────────────────────────────────────────

// Lyso-Gb3 derives from plasma Gb3 deacylation + tissue Gb3

dxdt_LGB3_PLM = kin_LGB3 * GB3_PLM - kel_LGB3 * LGB3_PLM;

// ── Inflammatory Index ────────────────────────────────────────────────────

// Inflammation driven by lyso-Gb3 (toxic signaling molecule)

dxdt_INFLAM = k_inflam_in * LGB3_PLM / 5.0

- k_inflam_out * INFLAM;

// ── eGFR (kidney function) ────────────────────────────────────────────────

// eGFR declines with persistent kidney Gb3 and proteinuria

double eGFR_loss_rate = k_eGFR_Gb3 * GB3_KID * (UPCR / 300.0 + 0.2);

dxdt_eGFR = -eGFR_loss_rate * eGFR;

// ── UPCR ─────────────────────────────────────────────────────────────────

// Proteinuria driven by podocyte Gb3 damage

double UPCR_driven = 40.0 * (GB3_KID / Gb3k_ss0);

dxdt_UPCR = 0.05 * (UPCR_driven - UPCR);

// ── LVMi ─────────────────────────────────────────────────────────────────

// Cardiac mass driven by Gb3 accumulation, constrained by minimum

double LVMi_driven = LVMi_min + k_LVMi_Gb3h * (GB3_HRT / Gb3h_ss0) *

(LVMi0 - LVMi_min);

double LVMi_capped = LVMi_driven < LVMi_min ? LVMi_min : LVMi_driven;

dxdt_LVMi = 0.02 * (LVMi_capped - LVMi);

// ── Neuropathic Pain ──────────────────────────────────────────────────────

// Pain depends on small fiber neuropathy driven by DRG Gb3

double E_eff = E_GalA / (E_GalA + 1.0); // normalized enzyme effect (0–1)

double Pain_target = Pain0 * (1.0 - 0.6 * E_eff) + 1.0; // floor = 1.0

dxdt_PAIN = 0.08 * (Pain_target - PAIN);

$TABLE

capture C_AGAB_ngmL = A_AGAB_C / (Vc_AGAB * BW);

capture C_AGAA_ngmL = A_AGAA_C / (Vc_AGAA * BW);

capture C_MIG_ugmL = C_MIG;

capture C_LUC_ugmL = C_LUC;

capture LYS_ERT_total = A_AGAB_LYS + A_AGAA_LYS;

capture Enzyme_activity_total = E_GalA;

capture GCS_inhibition_pct = GCS_inhib * 100;

capture Plasma_Gb3 = GB3_PLM;

capture Kidney_Gb3 = GB3_KID;

capture Cardiac_Gb3 = GB3_HRT;

capture LysoGb3 = LGB3_PLM;

capture Inflammation = INFLAM;

capture eGFR_val = eGFR;

capture UPCR_val = UPCR;

capture LVMi_val = LVMi;

capture BPI_Pain = PAIN;

# Compile the model

mod <- mcode("FabryDisease_QSP", fbr_code)

cat("Model compiled successfully.\n")

cat("Compartments:", length(mod@cmts), "\n")

print(mod)

##############################################################################

# TREATMENT SCENARIO DEFINITIONS

# Based on key clinical trials

##############################################################################

# S1: Natural History (Classic Fabry, no treatment)

# Reference: Schiffmann 2001, Mehta 2009 Eur J Clin Invest

# S2: Agalsidase Beta 1 mg/kg IV Q2W (Fabrazyme)

# Reference: Eng 2001 NEJM (FABRY-001); Banikazemi 2007 Ann Intern Med

# S3: Agalsidase Alfa 0.2 mg/kg IV Q2W (Replagal)

# Reference: Schiffmann 2001 Ann Intern Med; Schiffmann 2006 JAMA

# S4: Migalastat 150 mg QOD (oral, amenable mutations ~40% of patients)

# Reference: Germain 2016 NEJM (ATTRACT); Hughes 2017 Lancet

# S5: Pegunigalsidase Alfa 1 mg/kg IV Q4W (Elfabrio)

# Reference: Schiffmann 2021 JAMA (BRIGHT); Tøndel 2024

# S6: Lucerastat 1000 mg TID (SRT, GCS inhibitor, combination with ERT)

# Reference: Lenders 2022 Lancet DE (MODIFY trial)

# Simulation period: 5 years (1825 days)

SIM_END <- 5 * 365 # days

BW_KG <- 70

# Helper: Q2W dosing events (every 14 days)

make_q2w_events <- function(dose_ug, CMT_name, CMT_num, n_doses = 130) {

data.frame(

time = seq(0, by = 14, length.out = n_doses),

amt = dose_ug,

cmt = CMT_num,

evid = 1,

rate = dose_ug / 0.5 # 30-min infusion

)

}

# Helper: TID oral dosing

make_tid_events <- function(dose_ug, CMT_num, n_days = 365*5) {

times <- c(outer(0:(n_days-1) * 24, c(0, 8, 16), "+"))

data.frame(time = sort(times) / 24, # convert h->days? Use hours below

amt = dose_ug,

cmt = CMT_num,

evid = 1)

}

##############################################################################

# SCENARIO 1: Natural History (Classic Fabry, no treatment)

##############################################################################

cat("\n--- Scenario 1: Natural History ---\n")

s1_data <- data.frame(

ID = 1, time = 0, amt = 0, cmt = 1, evid = 0

)

s1 <- mod %>%

param(IS_CLASSIC = 1, IS_AMENABLE = 0) %>%

mrgsim(data = s1_data,

end = SIM_END, delta = 1,

carry_out = "evid") %>%

as.data.frame() %>%

mutate(Scenario = "S1: 자연경과 (치료 없음)", Year = time/365)

##############################################################################

# SCENARIO 2: Agalsidase Beta 1 mg/kg IV Q2W (Fabrazyme)

##############################################################################

cat("--- Scenario 2: Agalsidase Beta (Fabrazyme) 1 mg/kg Q2W ---\n")

# 1 mg/kg × 70 kg = 70 mg = 70,000 μg; Q2W

s2_ev <- ev(

amt = 70000 * 1000, # ng (model units)

cmt = 1, # A_AGAB_C

ii = 14, # every 14 days

addl = 129, # 130 doses total

rate = 70000 * 1000 / 0.5, # 30-min infusion (days)

evid = 1

)

s2 <- mod %>%

param(IS_CLASSIC = 1, IS_AMENABLE = 0, BW = 70) %>%

mrgsim(events = s2_ev,

end = SIM_END, delta = 1) %>%

as.data.frame() %>%

mutate(Scenario = "S2: 아갈시다제 베타 1mg/kg Q2W", Year = time/365)

##############################################################################

# SCENARIO 3: Agalsidase Alfa 0.2 mg/kg IV Q2W (Replagal)

##############################################################################

cat("--- Scenario 3: Agalsidase Alfa (Replagal) 0.2 mg/kg Q2W ---\n")

# 0.2 mg/kg × 70 kg = 14 mg = 14,000 μg; Q2W

s3_ev <- ev(

amt = 14000 * 1000, # ng

cmt = 4, # A_AGAA_C

ii = 14,

addl = 129,

rate = 14000 * 1000 / 0.5,

evid = 1

)

s3 <- mod %>%

param(IS_CLASSIC = 1, IS_AMENABLE = 0, BW = 70) %>%

mrgsim(events = s3_ev,

end = SIM_END, delta = 1) %>%

as.data.frame() %>%

mutate(Scenario = "S3: 아갈시다제 알파 0.2mg/kg Q2W", Year = time/365)

##############################################################################

# SCENARIO 4: Migalastat 150 mg QOD (Oral Chaperone, Amenable Mutation)

##############################################################################

cat("--- Scenario 4: Migalastat 150 mg QOD (amenable mutation) ---\n")

# 150 mg = 150,000 μg every 48 h

# Use IS_AMENABLE = 1 to activate pharmacological chaperone effect

s4_ev <- ev(

amt = 150000, # μg

cmt = 7, # A_MIG_GUT

ii = 2, # every 2 days (QOD)

addl = 911, # ~5 years

evid = 1

)

s4 <- mod %>%

param(IS_CLASSIC = 0, IS_AMENABLE = 1, BW = 70) %>% # late-onset + amenable

mrgsim(events = s4_ev,

end = SIM_END, delta = 1) %>%

as.data.frame() %>%

mutate(Scenario = "S4: 미갈라스타트 150mg QOD (샤페론)", Year = time/365)

##############################################################################

# SCENARIO 5: Pegunigalsidase Alfa 1 mg/kg IV Q4W (Elfabrio)

##############################################################################

cat("--- Scenario 5: Pegunigalsidase Alfa 1 mg/kg Q4W ---\n")

# 1 mg/kg × 70 kg = 70 mg; Q4W (every 28 days)

# PEGylation → prolonged t1/2; model as agalsidase beta with reduced CL

s5_ev <- ev(

amt = 70000 * 1000, # ng

cmt = 1,

ii = 28, # Q4W

addl = 64, # 65 doses

rate = 70000 * 1000 / 3.5, # 3.5-h infusion

evid = 1

)

s5 <- mod %>%

param(IS_CLASSIC = 1, IS_AMENABLE = 0, BW = 70,

CL_AGAB = 0.05, # lower CL due to PEGylation → t1/2 ~80h

Km_ERT = 0.8, # improved lysosomal targeting

Emax_ERT = 80.0) %>%

mrgsim(events = s5_ev,

end = SIM_END, delta = 1) %>%

as.data.frame() %>%

mutate(Scenario = "S5: 페구니알시다제 알파 1mg/kg Q4W", Year = time/365)

##############################################################################

# SCENARIO 6: ERT + Lucerastat Combination (SRT add-on)

##############################################################################

cat("--- Scenario 6: Agalsidase Beta + Lucerastat 1000 mg TID ---\n")

# ERT: same as S2 Q2W

# Lucerastat: 1000 mg TID; model as QD equiv (simplified for mrgsolve)

s6_ert <- ev(

amt = 70000 * 1000, cmt = 1,

ii = 14, addl = 129,

rate = 70000 * 1000 / 0.5, evid = 1

)

s6_luc <- ev(

amt = 1000000, # 1000 mg = 1,000,000 μg

cmt = 9, # A_LUC_GUT

ii = 1.0/3.0, # every 8h (TID)

addl = 5*365*3-1,

evid = 1

)

s6_ev <- c(s6_ert, s6_luc)

s6 <- mod %>%

param(IS_CLASSIC = 1, IS_AMENABLE = 0, BW = 70) %>%

mrgsim(events = s6_ev,

end = SIM_END, delta = 1) %>%

as.data.frame() %>%

mutate(Scenario = "S6: ERT + 루세라스탓 1000mg TID 병용", Year = time/365)

##############################################################################

# COMBINE ALL SCENARIOS

##############################################################################

all_scenarios <- bind_rows(s1, s2, s3, s4, s5, s6) %>%

filter(time >= 0)

cat("\n=== Simulation Complete ===\n")

cat("Scenarios:", length(unique(all_scenarios$Scenario)), "\n")

cat("Total rows:", nrow(all_scenarios), "\n")

##############################################################################

# VISUALIZATION

##############################################################################

theme_fabry <- theme_bw(base_size = 12) +

theme(

legend.position = "bottom",

legend.title = element_blank(),

strip.background = element_rect(fill = "#E3F2FD"),

plot.title = element_text(face = "bold", size = 14),

panel.grid.minor = element_blank()

)

scenario_colors <- c(

"S1: 자연경과 (치료 없음)" = "#616161",

"S2: 아갈시다제 베타 1mg/kg Q2W" = "#1565C0",

"S3: 아갈시다제 알파 0.2mg/kg Q2W" = "#0288D1",

"S4: 미갈라스타트 150mg QOD (샤페론)" = "#2E7D32",

"S5: 페구니알시다제 알파 1mg/kg Q4W" = "#7B1FA2",

"S6: ERT + 루세라스탓 1000mg TID 병용" = "#E65100"

)

# ── Plot 1: Plasma Lyso-Gb3 (most sensitive biomarker) ──────────────────

p1 <- all_scenarios %>%

ggplot(aes(Year, LysoGb3, color = Scenario)) +

geom_line(size = 1.2) +

scale_color_manual(values = scenario_colors) +

labs(title = "혈장 Lyso-Gb3 (파브리병 핵심 바이오마커)",

x = "연도 (Years)", y = "Lyso-Gb3 (μg/L)") +

geom_hline(yintercept = 2.0, linetype = "dashed", color = "gray50") +

annotate("text", x = 4.5, y = 2.5, label = "정상 상한값 (<2 μg/L)",

size = 3.5, color = "gray50") +

theme_fabry

print(p1)

# ── Plot 2: eGFR over time ──────────────────────────────────────────────

p2 <- all_scenarios %>%

ggplot(aes(Year, eGFR_val, color = Scenario)) +

geom_line(size = 1.2) +

scale_color_manual(values = scenario_colors) +

labs(title = "eGFR 경과 (신장 기능 보존 효과)",

x = "연도 (Years)", y = "eGFR (mL/min/1.73m²)") +

geom_hline(yintercept = 60, linetype = "dashed", color = "red") +

annotate("text", x = 4.5, y = 62, label = "CKD G3 경계 (60)",

size = 3.5, color = "red") +

theme_fabry

print(p2)

# ── Plot 3: LVMi ──────────────────────────────────────────────────────

p3 <- all_scenarios %>%

ggplot(aes(Year, LVMi_val, color = Scenario)) +

geom_line(size = 1.2) +

scale_color_manual(values = scenario_colors) +

labs(title = "LVMi — 좌심실 질량 지수 (심장 비후 역전)",

x = "연도 (Years)", y = "LVMi (g/m²)") +

geom_hline(yintercept = 115, linetype = "dashed", color = "blue") +

annotate("text", x = 4.5, y = 118, label = "남성 정상 상한값 (115 g/m²)",

size = 3.5, color = "blue") +

theme_fabry

print(p3)

# ── Plot 4: Neuropathic Pain ──────────────────────────────────────────

p4 <- all_scenarios %>%

ggplot(aes(Year, BPI_Pain, color = Scenario)) +

geom_line(size = 1.2) +

scale_color_manual(values = scenario_colors) +

ylim(0, 10) +

labs(title = "신경병성 통증 (BPI-SF, 0–10)",

x = "연도 (Years)", y = "통증 점수 (0–10)") +

theme_fabry

print(p4)

# ── Plot 5: Enzyme Activity ───────────────────────────────────────────

p5 <- all_scenarios %>%

filter(Year <= 0.5) %>% # short term to see PK peaks

ggplot(aes(Year * 365, Enzyme_activity_total, color = Scenario)) +

geom_line(size = 1.1) +

scale_color_manual(values = scenario_colors) +

labs(title = "α-갈락토시다제 A 효소 활성 (초기 6개월)",

x = "일수 (Days)", y = "α-Gal A (nmol/h/mg)") +

theme_fabry

print(p5)

# ── Plot 6: Multi-panel summary ───────────────────────────────────────

library(patchwork)

p_combined <- (p1 + p2) / (p3 + p4) +

plot_annotation(

title = "파브리병 QSP 모델: 치료 시나리오 비교 (5년 시뮬레이션)",

subtitle = "ERT vs 미갈라스타트 샤페론 vs SRT 병용 — 신장·심장·신경계 복합 엔드포인트",

theme = theme(plot.title = element_text(face = "bold", size = 16))

)

print(p_combined)

##############################################################################

# SUMMARY TABLE: 5-YEAR OUTCOMES

##############################################################################

summary_table <- all_scenarios %>%

filter(abs(Year - 5) < 0.01) %>%

group_by(Scenario) %>%

summarise(

eGFR_5yr = round(mean(eGFR_val), 1),

eGFR_change = round(mean(eGFR_val) - 90.0, 1),

LVMi_5yr = round(mean(LVMi_val), 1),

LVMi_change = round(mean(LVMi_val) - 148.0, 1),

LysoGb3_5yr = round(mean(LysoGb3), 2),

Pain_5yr = round(mean(BPI_Pain), 1),

.groups = "drop"

)

cat("\n=== 5-Year Outcome Summary ===\n")

print(summary_table)

cat("\n=== Parameter Calibration Notes ===\n")

cat("ERT PK: Vedder AC (2007) J Am Soc Nephrol (agalsidase beta)\n")

cat(" Schiffmann R (2001) Ann Intern Med (agalsidase alfa)\n")

cat("Migalastat PD: Germain DP (2016) NEJM (ATTRACT), EC50 ~0.2 μg/mL\n")

cat("Pegunigalsidase: Schiffmann R (2021) JAMA (BRIGHT), extended t1/2 ~80h\n")

cat("Lucerastat: Lenders M (2022) Lancet DE (MODIFY), GCS IC50 ~0.18 μg/mL\n")

cat("eGFR natural history: Warnock DG (2012) Am J Kidney Dis,\n")

cat(" -3 to -12 mL/min/1.73m²/yr without treatment\n")

cat("LVMi: Weidemann F (2009) Circulation, ERT slows/reverses hypertrophy\n")

cat("Lyso-Gb3: Aerts JM (2008) PNAS, most sensitive biomarker for response\n")

FilesExpand file tree

fbr_mrgsolve_model.R

Latest commit

History

fbr_mrgsolve_model.R

File metadata and controls