A quantitative systems pharmacology model of traumatic spinal cord injury, built around a single question: why has almost every neuroprotective drug that worked in the laboratory failed in the clinic? The model's answer is not "the drugs don't work". It is that the cascade they target is over in hours, that the endpoint they are measured on is a steep saturating function of the tissue they save, and that the patients easiest to enrol are exactly the patients in whom a real tissue benefit is arithmetically undetectable.

Scale: 238-node / 18-cluster / 378-edge mechanistic map · 43-ODE mrgsolve model with 160 annotated parameters · 18 prebuilt scenarios · 12-tab Shiny dashboard · 70 references (PMID-verified where the record could be confirmed, PubMed search links otherwise, with an explicit list of what is not sourced).

Every quantitative claim in this README is produced by python3 sci_reference_model.py, a dependency-free re-implementation of the same 43-state system, and the full run is saved in sci_reference_output.txt. No number here is asserted from memory.

Traumatic spinal cord injury begins with an irreversible mechanical insult — contusion, persistent compression, axonal shear — that no drug will ever reach, and continues as a secondary injury cascade that is, in principle, entirely addressable. Residual compression plus cord swelling raise intraspinal pressure (ISP) while neurogenic shock lowers mean arterial pressure, and the resulting collapse in spinal cord perfusion pressure (SCPP = MAP − ISP) drives ischemia → ATP failure → SUR1-TRPM4/AQP4-mediated cytotoxic edema, which raises ISP further: a compression–edema–ischemia vicious cycle. Inside it, glutamate floods the extracellular space within minutes, Ca²⁺ overload activates calpain and opens the mitochondrial permeability transition pore, peroxynitrite and Fenton chemistry peroxidize membrane lipids, neutrophils arrive at ~24 h and macrophages over days to weeks, and caspase-3 kills neurons and oligodendrocytes for a fortnight after the accident. What is left is a cystic cavity walled off by a CSPG-rich glial scar, a demyelinated and thinned rim of spared white matter, and a spinal cord below the lesion that slowly reorganizes — restoring reflexes that produce spasticity, detrusor overactivity and autonomic dysreflexia while a closing critical period allows rehabilitation to recruit whatever descending drive survived. Function at one year is not a property of the injury alone: it is the product of what survived, how well it conducts, and how much of it the patient's circuits learned to use.

1. Primary mechanical injury (t = 0, irreversible, not a drug target): mechanism, impact energy, compression, contusion, shear, neurological level, initial AIS grade, central hemorrhagic necrosis, primary axonal loss, spared subpial rim
2. Systemic haemodynamics & cord perfusion: neurogenic shock, sympathectomy, bradycardia, MAP, vasopressor support, the MAP ≥ 85 mmHg protocol, ISP, SCPP, autoregulation loss, SCBF, ischemic penumbra, reperfusion, duraplasty
3. Cord swelling & ionic dysregulation: ATP depletion, Na⁺/K⁺-ATPase failure, SUR1–TRPM4, AQP4, cytotoxic and vasogenic edema, BSCB breakdown, MMP-9, tight-junction loss, rostral ascending edema
4. Glutamatergic excitotoxicity: vesicular release, EAAT/GLT-1 reversal, NMDA/AMPA-kainate, persistent Na⁺ current, Ca²⁺ overload, calpain, α-II-spectrin breakdown, mPTP, cytochrome-c, nNOS
5. Oxidative stress: superoxide, peroxynitrite, hydroxyl radical, free iron, lipid peroxidation, 4-HNE/MDA, protein nitration, ferroptosis, Nrf2/HO-1, SOD/catalase
6. Neuroinflammation: DAMPs/HMGB1, TLR4–NF-κB, NLRP3, IL-1β/TNF-α/IL-6, CXCL1/CCL2, neutrophils, monocytes, microglia, the M1/M2 balance, lymphocytes, debris clearance, SCI-induced immune depression
7. Cell death programs: caspase-8/9/3, Bax/Bcl-2, necroptosis, autophagy, neuron death, oligodendrocyte apoptosis, Wallerian degeneration, retrograde dieback, cystic cavitation, late syringomyelia
8. Demyelination & remyelination: myelin loss, conduction block, juxtaparanodal K⁺ exposure, Nav redistribution, OPC/NG2 proliferation, Schwann-cell remyelination, thin sheaths
9. Glial scar & regeneration failure: reactive astrogliosis, fibrotic core, CSPG, PTPσ/LAR, Nogo-A, MAG/OMgp, NgR1/p75/LINGO-1, RhoA–ROCK, growth-cone collapse, PTEN/mTOR, KLF/RAG silencing, perineuronal nets
10. Spared tracts, plasticity & circuit reorganization: CST and reticulospinal sparing, propriospinal relays, effective descending drive, collateral sprouting, synaptogenesis, the critical period, BDNF/NT-3, locomotor CPG, cortical remapping, maladaptive plasticity
11. Methylprednisolone PK/PD including start time, cumulative exposure and the complication set (pneumonia, GI haemorrhage, hyperglycaemia, myopathy)
12. Riluzole PK/PD (CYP1A2, persistent Na⁺-current block, release inhibition, enhanced uptake, ALT elevation)
13. Minocycline & glibenclamide PK/PD (M1 blockade, MMP inhibition, caspase inhibition; SUR1–TRPM4 blockade, hypoglycaemia)
14. Regenerative & neuromodulation interventions: intrathecal anti-Nogo-A, Rho-ROCK inhibition, chondroitinase, activity-based rehabilitation, body-weight-supported and robotic gait training, epidural stimulation, FES cycling, cell therapy, decompression timing
15. Symptomatic pharmacology: baclofen (and its weakness trade-off), pregabalin, antimuscarinics/onabotulinumtoxinA, α-blockade
16. Secondary systemic complications: reflex reorganization, spasticity, neuropathic pain, neurogenic bladder and DSD, autonomic dysreflexia, respiratory failure, sublesional osteoporosis, pressure injury, VTE, cardiometabolic decline, heterotopic ossification
17. Clinical endpoints: ISNCSCI motor and sensory scores, AIS conversion, SCIM III, WISCI II, 10-m walk, GRASSP, pain NRS, Ashworth, MRI lesion volume and tissue bridges, MEP/SSEP, CSF NfL/GFAP, QoL, mortality/LOS
18. Quantitative concepts encoded in the ODEs, drawn as explicit notes on the map: the cascade flux Φ(t), the window integral, the threshold mapping, surrogate dissociation, the two clocks, and recovery-creates-pathology

Two structural choices carry most of the model's behaviour:

• AXON only ever falls; PLAST is what rises. Effective descending drive is CONN = AXON × myelin efficiency × conduction factor × (1 + plasticity + stimulation gain). Anatomy and function are separate variables, which is the only way a plasticity-directed therapy can change outcome without touching the imaging surrogate — and vice versa.
• Only supra-threshold death signalling destroys tissue (APOP_THR, ISCH_THR). Without this the system drains axons forever at its own tiny chronic steady state; with it, the chronic state is genuinely quiescent and late decline has to be put in deliberately.

Standard care in the model = decompression at 24 h, rehabilitation intensity 0.5 from day 14. For a C5 lesion, initial → 1-year ISNCSCI total motor score:

The day-3 nadir is not tissue loss — it is reversible conduction block from edema and ischemia. That is why patients improve from their admission exam at all, and why the admission exam systematically understates the surviving anatomy. About 90% of the modelled change is complete by 3–6 months.

The NASCIS methylprednisolone regimen (30 mg/kg bolus + 5.4 mg/kg/h × 23 h) given at increasing delays, in a C5 AIS C patient:

Benefit falls 3.1-fold from 1 h to 24 h. The complication index is identical to three decimal places at every delay, because efficacy is ∫ E(t)·Φ(t) dt against a cascade flux that decays, while exposure is ∫ C(t) dt against nothing. Doubling the infusion to 48 h buys 1.39× the complication index for +0.42 motor points — the drug has already spent its window by the time the second day of infusion begins. This is the entire NASCIS story reproduced from mechanism rather than from post-hoc subgroups: in the same patient at AIS A, the benefit is +0.31 at 1 h and +0.08 at 24 h — a real drug effect that no trial could measure.

MOTOR is a steep saturating Hill function of CONN (TH50 = 0.22, HILL = 2.2). The value of doubling a patient's spared drive:

Run the same drug across grades and the dissociation becomes explicit. Riluzole (14 d, RISCIS regimen) reduces the lesion index by 8.7–9.4% in every grade — essentially grade-independent, as a tissue-level mechanism should be — while its effect on the endpoint is not:

A neuroprotectant can therefore be mechanistically successful and clinically negative at the same time, and the difference between those two outcomes is who was enrolled, not whether the drug works. AIS A patients are the easiest to recruit and consent, have the least ambiguous baseline, and sit exactly where the derivative of the endpoint with respect to the mechanism is smallest.

The neuroprotection clock runs in hours. The plasticity clock runs in weeks and is shut by the scar (CRIT decays at K_CRIT·(1 + K_CRIT_SCAR·CSPG)). The same rehabilitation dose, moved 106 days later (C5 AIS C):

+8.2 motor points for the same intervention delivered earlier — a window effect that has nothing to do with the acute cascade and cannot be recovered by doing more rehabilitation later.

Meanwhile the sublesional reflex reorganization that underlies functional recovery is itself the source of the chronic complications. The same bladder distension applied at different times after a T4 injury:

A harmless nursing event in week 1 is a hypertensive emergency at 6 months, and nothing about the lesion changed in between. The same variable drives Ashworth from 0.0 → 2.7 and pain NRS from 0.2 → 3.6 over the first three months.

C5 AIS A, 1 year, versus standard care (lesion index 0.323, motor 12.1):

Two things worth reading twice. Haemodynamic management, which involves no new molecule, dominates every drug arm in this model, because it acts on the SCPP term that gates the entire cascade for as long as the cord is swollen — that is a testable prediction, and SCPP/ISP are captured as outputs so it can be argued with. And symptomatic pharmacology is not free: baclofen lowers Ashworth through the same GABA-B mechanism that produces weakness, so the modelled motor score falls by 2.6 points. The trade-off is structural, not a side effect.

For an incomplete injury, epidural stimulation from day 120 is the largest single effect anywhere in the model: AIS C 51.3 → 61.9 motor points, with no change whatsoever in lesion index or axon sparing — it amplifies residual drive rather than saving tissue.

Sublesional bone loss (−22% at one year without loading) is driven by absent weight-bearing, not by the lesion — so it is the one endpoint that FES and standing improve without any change in descending drive. It is also the reason the model tracks BMD at all: an intervention that never moves the motor score can still move a hard outcome.

# mrgsolve
library(mrgsolve)
mod <- mread_cache("sci_mrgsolve_model.R")

# all 18 scenarios (helpers are documented at the bottom of the model file)
df <- SCI_simulate_scenarios(365)

# the two window experiments
SCI_mp_window(ais = 3)            # methylprednisolone start time, AIS C
SCI_decompression_sweep(ais = 3)  # time to surgical decompression, AIS C
SCI_ad_challenge()                # autonomic dysreflexia, T4, by trigger day

# Shiny dashboard (12 tabs; model file must be in the same directory)
shiny::runApp("sci_shiny_app.R")

# dependency-free reference implementation — reproduces every number above
python3 sci_reference_model.py          # full report (~3 min)
python3 sci_reference_model.py --quick  # scenarios only

Render the map:

dot -Tsvg sci_qsp_model.dot -o sci_qsp_model.svg
dot -Tpng -Gdpi=150 sci_qsp_model.dot -o sci_qsp_model.png

Patient & injury profile · Drug PK · Secondary cascade · Neuroinflammation · Tissue & imaging surrogates · Neurological recovery · Therapeutic-window explorer · Surrogate-vs-endpoint explorer · Secondary complications · Scenario comparison · Steroid safety ledger · References.

The two bolded tabs carry the argument: the first sweeps intervention delay and plots benefit against the complication index on the same axes; the second marks where the current patient sits on the Hill curve and shows what a doubling of their spared drive would be worth in ISNCSCI points.

Read sci_references.md § 16 before reusing anything here. In short: most rate constants are shape-matched to described time courses rather than fitted; TH50/HILL are a calibrated modelling choice, so the direction of the surrogate–endpoint argument is robust but the specific point estimates are not; the lesion index is a bounded 0–1 index with a nominal mL conversion and must not be compared numerically to MRI volumetry; drug PK is dose-proportional rather than physiological; the complication index orders regimens by exposure and is not a probability; LEVEL_IDX is a lookup, not a segmental cord, so central-cord and Brown-Séquard syndromes cannot be represented; and syringomyelia, pressure injury, VTE, heterotopic ossification and the cardiometabolic trajectory appear on the map but have no ODEs.

Research and education only. Not a substitute for clinical judgment, and not suitable for regulatory use.