Getting Started

A CT scan shows an organ's large vessels and nothing smaller. Vasculature.jl grows the rest of the tree, down to the terminal arterioles, then lets physics decide how blood moves through it. It sends a contrast injection through the whole body and into that tree, and writes a CT phantom in which the true perfusion of every voxel is known.

This page takes you from installing the package to your first phantom. It shows how the heart, the lungs and the leg were built, and where to start for a new organ or a new disease. It also lists what went wrong for us along the way, so it does not go wrong for you.

Why Grow a Tree?

CT perfusion measures blood flow from how fast contrast fills the tissue. To test a method like that, you need tissue whose true flow you know, and patients never tell you. Digital phantoms such as XCAT draw an organ's anatomy beautifully, but they stop at the vessels a scanner resolves: in the heart, the coronary arteries and veins on its surface.

So we grow the rest. From XCAT's coronary arteries the heart grows 957,045 terminals, one for every small region of myocardium a terminal arteriole feeds. Flow then follows from first principles. Nothing in the pipeline is told what the perfusion should be: in the healthy heart it comes out at 0.86 mL/min/g at rest and 3.85 under adenosine, a flow reserve of 4.5. A CT simulator such as BasisSimulator.jl scans the phantom, and every measurement can be checked against the truth, voxel by voxel.

One Code Path for Every Organ

Four functions carry an organ from its imaged vessels to a phantom:

Each stage takes a spec: GrowthSpec, FlowSpec, TransportSpec, PhantomSpec. A spec is a plain Julia struct, and every value in it is a measurement with a source. The heart is HEART_GROWTH, HEART_FLOW, HEART_TRANSPORT and HEART_PHANTOM; the lungs are the same four structs with the lungs' numbers. A variant is one keyword: FlowSpec(HEART_FLOW; p_in = 70).

Nothing in the package asks which organ it is running. A brain will not be a new code path. It will be four new specs, and a notebook that prepares the brain's geometry. This is the most important design decision in the package, and the one to protect when you add to it.

Install

Julia 1.12 or later. Growth runs on the CPU or on any GPU.

using Pkg
Pkg.add(url = "https://github.com/MolloiLab/Vasculature.jl")
Pkg.add("Metal") # Apple silicon; CUDA for NVIDIA, AMDGPU for AMD, or nothing for the CPU

The XCAT phantom is licensed. Its files, and everything made from them (meshes, label volumes, grown trees, phantoms, scans), stay on lab storage. They never go into a public repository, an issue or a chat.

How We Work: Notebooks in SpaceStation

Every organ is built in a Pluto notebook in docs/notebooks/: 00_cardiac_pipeline, 01_pulmonary_pipeline, 02_limb_pipeline. A notebook is a plain Julia file. It prepares the organ's geometry, calls the four stages, checks the results and draws them. The same file is the documentation: the Examples pages are these notebooks.

We run them in SpaceStation, a workspace for Pluto notebooks that people and agents share:

spacestation --agents-md docs/notebooks

It starts in lazy mode. Nothing runs until you ask, and when a file changes, the cells it touches turn stale instead of running. Outputs are cached beside the notebook and survive restarts. That matters here, because a growth can take minutes and a phantom ten: you run what changed, and nothing else.

To work with an agent, start Claude Code in a terminal inside SpaceStation. It edits the notebook as a file, you watch the changed cells go stale, and it runs only those with pluto-collab run <notebook>.jl --stale. It never reruns the whole notebook, and pluto-collab restart is only for a worker that died.

Every heavy step in a notebook goes through cached(name) do … end: the first run computes and saves, and every run after loads in seconds. A notebook includes like any Julia file, so a script on a cluster can include it and reuse every function and cache in it.

Where the Flow Comes From

The rule that shapes everything else is that flow is never prescribed. Growth knows nothing about flow targets. The flow model knows two pressures, the trees' geometry and the physiology of a vessel wall, and the ground truth is whatever comes out. Here is how the heart gets its flow.

1. Maximal dilation is calibrated once. Under adenosine every arteriole is fully relaxed, so resistance alone sets the flow. Grown arterioles are longer and thinner than real ones (growth runs each terminal to the middle of its territory, where real tissue has a capillary network), so their raw resistance is too high. One factor, κ, scales the microvessels so that the whole bed's resistance matches healthy people measured by PET (r_total_dilated). Their resistances relative to each other, and so the spread of flow from voxel to voxel, still come from the tree.
2. Rest follows the oxygen balance. At rest the heart sets its flow to its oxygen need. Every arteriole below 150 µm is a wall in force balance. Its tone follows wall tension (the myogenic response) and the oxygen state of the beds it feeds (ATP released by desaturating red cells, conducted upstream). One tone constant is set so the blood leaving the heart has the measured oxygen extraction, 0.73. After that, each arteriole finds its own diameter.
3. Everything else is a prediction. The distribution of flow, the pressure–flow curve, transit times and the response to a stenosis are not fitted to anything. They are what the physics gives, and the tests check them against experiments.

Oxygen enters in two places and only two. The tissue's O₂ demand (mvo2) and extraction (oef) set resting flow through Fick's principle. Fick also sets how many terminals to grow: flow_ml_min, the organ's flow, over the flow of one terminal. The lungs use neither. They are not fed by the blood they carry, so their tone is passive, their terminal count is their structure (one per acinus), and their resistance is a prediction too.

Every value in a spec has its source in test/reference/, and the tests compare the model with measurements it was never fitted to. When the two disagree, we say so and leave it. The lungs carry 3.45 L/min, a quarter below the measured mean, because their grown veins are narrow; we report that rather than tune it away. A phantom tuned to look right cannot test anything.

Walkthrough: the Heart

The whole heart is the 00_cardiac_pipeline notebook. These are its four calls, with what each one does.

Grow

The tissue is XCAT's ventricular walls on a 1 mm grid. The inlets are XCAT's coronary arteries as trees of centrelines, and the chambers are blockers. Every imaged vessel first gets a short branch into the wall every 2.8 mm (a penetrating stub), so growth starts along its whole length and not only at its few imaged ends.

using Vasculature, Metal

spec = GrowthSpec(HEART_GROWTH; device = MtlArray)
tree = grow(spec, "heart.vol", ["lca.xml", "rca.xml"]; blockers = (chambers, origin_mm, 0.4))
veins = grow(GrowthSpec(HEART_VEIN_GROWTH; device = MtlArray), "heart.vol", ["veins.xml"])

HEART_GROWTH grows 957,045 terminals in about 40 seconds on a laptop GPU. Growth is HCT-CCO's staged constrained constructive optimization: each new terminal joins the tree where it adds the least blood volume, so the tree fills the tissue the way real trees do.

Flow

F = flow(HEART_FLOW, vtree(tree), vtree(veins), tissue) # tissue: 3 × n voxel centres (mm)
F.rest.resist.mbf, F.stress.resist.mbf # mL/min/g at rest and at maximal dilation

Transport

An injection into the right arm travels through a whole-body circulation (the venae cavae, the right heart, the lungs, the left heart and every organ bed), into the grown coronary tree, through every capillary bed and vein, and around again.

injection = Body.Injection(; phases = [(35, 5, 350), (35, 5, 105)]) # mL, mL/s, mg I/mL: contrast, then a 30:70 chaser
T = transport(HEART_TRANSPORT, F.rest, injection; times = HEART_PHANTOM.times)

Phantom

Every 0.2 mm voxel of XCAT's chest is classified: a tissue voxel, a vessel segment, a body pool (the aorta, a ventricle) or a static material. Each voxel also gets a compartment, the region ground truth is reported for.

phantom(HEART_PHANTOM, [("rest", F.rest, Trest), ("stress", F.stress, Tstress)], classify, grid;
    F.pairs, F.terr, nvox, out = "healthy", roi, roi_names, roi_groups,
    tracking = ["descending aorta" => (x, y, z)])

ph = DynamicPhantom.load("healthy"; rois = true)
P, materials = DynamicPhantom.bs_phantom(ph, BasisSimulator; state = "rest", t = 25)

The folder holds every voxel's label and compartment, each label's iodine every second, the ground-truth perfusion map, and an iodine ledger that must balance: what was delivered equals what is held plus what left through the veins. The build fails if it does not.

Three Organs, One Recipe

The heart, the lungs and the leg use the same four calls. What changes is the geometry each notebook prepares and the numbers in each spec.

The leg shows what happens when an organ is large. Murray's law says a tree's terminal radius is r_root · N^(−1/3), so reaching 10 µm from a 3.89 mm femoral artery takes 58.7 million terminals, too many for one growth. So the leg grows twice: first the whole leg down to 10 mL patches, then each patch down to 10 µm on its own. Every radius in both passes follows Murray's law from the terminals it feeds, so the result is one tree. The patches share no vessels, so they run in parallel.

Building a New Organ

Say you want a brain. Here is the order to work in, and where each piece lives.

1. Start the notebook. Copy the closest one (the heart, for an O₂-regulated organ) to docs/notebooks/03_brain_pipeline.jl, keep its sections, and open the folder in SpaceStation.
2. The tissue. A voxel mask of the tissue to perfuse, written with write_vol (or mesh_domain from a surface). For a brain, grey and white matter. A FlowSpec describes one tissue, so decide early whether they are one or two.
3. The inlets. The vessels the scan shows, as trees: XCAT.centerline turns a surface into a centreline, XCAT.base_forest or XCAT.tube_forest joins them, and XCAT.write_base_gxl writes them for grow. For a brain, both internal carotids and both vertebral arteries. Add penetrating stubs (XCAT.penetrating_stubs), and let the imaged lumens be crossed without taking terminals (XCAT.vessel_tube with NOTERM).
4. The blockers. Anything a vessel must not cross. For a brain, the ventricles and the skull.
5. The count. Fick for an organ its blood feeds (flow_ml_min: an older adult's 510 mL/min gives 3.6 million terminals), a structural count if the organ's units are known, or Murray from a target terminal size. If it is too large for one growth, grow in two passes, as the leg does.
6. The growth spec. GrowthSpec(HEART_GROWTH; tissue_ml, root_caps, gamma, flow_ml_min). Choose γ for the organ, then check it: the imaged vessels should come out near their imaged radii, and no capillary bed should end up with its venous pressure above its arterial one.
7. The flow spec. FlowSpec(HEART_FLOW; …) with the organ's own values: pressures, blood volume and its split, capillary size, O₂ demand and extraction, resistance at maximal dilation, permeability. Write each source into test/reference/. Then solve and check: no reversed beds, resting flow and venous saturation within the measured range, tone not saturated.
8. Transport. The body model has organ slots for the heart and the lungs. A brain, a liver or a leg first needs its own slot, with its share of the cardiac output taken from the grown tree's flow, as coronary_share and organ_flow do. That is package work: open an issue first.
9. The phantom. A classifier that says what every voxel is, compartments (each also grouped with the conduit vessels that feed and drain it), and a bolus-tracking point inside an imaged artery.
10. Into the package. When the notebook works, move its specs into src/organs.jl and src/perfusion/organs.jl, move any helper the next organ would need into src/, add tests, and add the notebook to the docs.

Building a New Disease

A disease is a change to a healthy organ, solved again with the healthy organ's tone constant so its arterioles respond on their own. The heart's 90% LAD stenosis is the example:

lesion = [(segment, area_fraction, length_mm) for …] # the first 10 mm of the LAD
stress = Flow.circuit(arteries, F.terr; state = :stress, veins, vleaf = F.vleaf, phys = HEART_FLOW,
    reference = F.stress, kappa = F.stress.resist.kappa, stenosis = lesion)
rest = Flow.circuit(arteries, F.terr; state = :rest, veins, vleaf = F.vleaf, phys = HEART_FLOW,
    reference = F.rest, kappa = F.stress.resist.kappa, stenosis = lesion, demand = F.demand)

The lesion is graded on XCAT's imaged lumen. Its pressure loss follows the viscous and separation terms of Young and Tsai, and the arterioles downstream dilate until they run out of reserve. The phantom then paints plaque (fibrous tissue, a lipid core, calcium) where the lesion narrows the lumen, and the scenario is saved as its own folder.

The lungs' occlusion is the same call with a tiny area fraction (10⁻³) on the artery that feeds most of one lobe. More stenoses, serial lesions, microvascular disease or a lower pressure are each a geometry change or a spec field on the same calls. Keep the healthy run as the reference, and give each scenario its own folder.

Rough Edges, and How We Got Past Them

Each of these cost us hours or days. They are written as what you will see, why it happens, and what to do. Open the ones that match what you are looking at.

Geometry and inputs

Growth

Flow

Transport and timing

Phantom and painting

Notebooks and workflow

The Rules We Keep

• One code path. One function per stage, the same for every organ. Organs differ only in their specs and their geometry; nothing in src/ branches on the organ.
• Specs are Julia structs. Presets live in code, every value cited, with its source in test/reference/. No configuration files.
• Flow is pressure-driven. Ground truth is what the physics gives. Measurements are comparators, never targets, and a disagreement is reported, not tuned away.
• Each organ has its own physiology. Never carry one organ's values, thresholds or exponent into another.
• Grow real terminals. Reach terminal-arteriole size with more terminals, never with a lower γ. Growth uses HCT-CCO's staged schedule only.
• A disease keeps the healthy tone. Solve it with the healthy circuit as its reference, so the arterioles respond on their own.
• Notebooks orchestrate; the package holds the logic. A helper a second organ needs moves into src/. Superseded code is deleted, not archived.
• Small numerical differences are fine. CPU, Metal and CUDA runs differ in the last digits. Do not chase bit-exact results.
• Licensed data stays on lab storage. XCAT, and everything made from it, never goes anywhere public.

Where to Go Next

• The Examples: the heart, the lungs and the leg, each as a notebook you can open in SpaceStation.
• The API: every exported function and spec, with its docstring.
• The source: src/spec.jl (growth), src/perfusion/flow.jl (flow), src/perfusion/transport.jl and src/perfusion/body.jl (transport), src/phantom.jl (phantom), src/xcat/ (XCAT inputs), src/organs.jl and src/perfusion/organs.jl (the presets).

For Your Agent

Paste this into any coding agent before it builds on the package. It holds the goals, the rules and the rough edges above in the form an agent needs.

You are working on Vasculature.jl (github.com/MolloiLab/Vasculature.jl), a Julia package that grows complete vascular trees from imaged vessels, solves blood flow from pressures, transports contrast through a whole-body model, and writes CT perfusion phantoms with known ground truth. Docs: molloilab.github.io/Vasculature.jl.

PIPELINE (one code path for every organ; organs differ only in specs and geometry)
- grow(spec::GrowthSpec, domain, inlets; blockers, territories) -> engine; vtree(engine) — src/spec.jl
- flow(spec::FlowSpec, arteries, veins, tissue_centres) -> (; rest, stress, terr, pairs, vleaf, demand); for a disease, Flow.circuit(...; reference, kappa, stenosis, demand) — src/perfusion/flow.jl
- transport(spec::TransportSpec, circuit, injection; times) — src/perfusion/transport.jl, body.jl
- phantom(spec::PhantomSpec, scans, classify, grid; pairs, terr, nvox, out, roi, roi_names, roi_groups, tracking) -> DynamicPhantom bundle — src/phantom.jl
- Presets: HEART_*, LUNG_* in src/organs.jl and src/perfusion/organs.jl; sources in test/reference/.
- XCAT inputs: src/xcat/ (centerline, base_forest, tube_forest, graft, drain, penetrating_stubs, vessel_tube, vessel_bridges, scale_tree, great_vessel_transit).
- Organ notebooks: docs/notebooks/00_cardiac_pipeline.jl, 01_pulmonary_pipeline.jl, 02_limb_pipeline.jl.

RULES
- Never add organ-specific branches to src/. A new organ = new GrowthSpec/FlowSpec/TransportSpec/PhantomSpec + a notebook that prepares its geometry. Helpers a second organ needs move into src/ with tests.
- Specs are Julia structs; every value measured and cited (test/reference/*.toml). No config files.
- Flow is never prescribed. Ground truth is what the physics gives. Do not tune growth or specs to hit a measurement; report disagreements with the reason (e.g. lung flow 3.45 vs 4.8 L/min).
- Each organ has its own physiology; never reuse another organ's values, thresholds or gamma.
- Reach real terminal size (~10 µm) with more terminals (two-pass growth if needed), never by lowering gamma. Growth uses the HCT-CCO staged schedule only.
- Diseases: re-solve with the healthy circuit as `reference` (keeps the tone constant); one folder per scenario.
- Never change a perfusion measurement method to agree with the phantom; fix the physics.
- Delete superseded code; no archives or .bak files. Float noise across CPU/Metal/CUDA is fine.
- XCAT data and anything derived from it is licensed: never commit, publish or paste it.

WORKFLOW (Pluto notebooks in SpaceStation, lazy mode)
- Notebooks are plain .jl files; edit them as files. Keep a blank line before every "# ╔═╡" cell marker (Pluto drops a cell's last character otherwise). Update the "# ╔═╡ Cell order:" block when adding cells; fold table/figure/markdown cells (╟─), keep TableOfContents() shown.
- Run only what changed: `pluto-collab run <nb>.jl --stale`; check `pluto-collab status <nb>.jl`. Never rerun a whole notebook; `pluto-collab restart` only for a dead worker.
- Wrap heavy steps in cached(name) do ... end; delete a cache after changing a struct it holds.
- Put loops over notebook globals in functions; draw many segments with one linesegments!.
- Docs: render with docs/extract_all.jl, verify with docs/verify_notebook_exports.jl, and commit the notebook and its export together (a live session can re-save the notebook after rendering).

ROUGH EDGES ALREADY SOLVED (do not rediscover)
- Growth only from imaged ends -> XCAT.penetrating_stubs (2.8 mm x size scale).
- Branches collide with imaged pieces -> resample centrelines every 2.8 mm x size scale.
- Imaged vessels outside the tissue -> vessel_tube / vessel_bridges as NOTERM voxels.
- XCAT anatomy off-normal (aorta ~3 SD narrow) -> time positional vessels from the body model (scale_tree); check every vessel against measurements.
- Beds reverse with the coronary gamma on a large organ -> gamma per organ (leg: 3).
- Terminals too large -> more terminals; r_terminal = r_root * N^(-1/gamma).
- Two-pass patches too thin / seams -> Murray radii from final counts (r = 10 µm * n^(1/3)); 3 mm NOTERM margin into neighbouring patches.
- NaN pressures on a new organ -> the organ's own FlowSpec.
- Stress perfusion low -> PS rises with flow (Renkin-Crone, PS = 0.22 + 0.35 MBF), fb_dilated.
- Body vs organ flow -> coronary_share (heart), organ_flow (lungs); new organs need a body slot (ORGAN_SLOTS is (:coronary, :lung) today) with the same rule.
- AIF timing -> positional great-vessel transit; ideal compartment includes conduit vessels.
- Tracking curve wrong -> tracked segment has its own class at blood concentration; AIF ROI inside it.
- Ledger must balance within 1%; a failure is a classifier/compartment bug.
- XCAT labels: ids reused by distant organs (restrict by region); label 7 = airway wall, 11 = air; the main pulmonary artery is painted as RV blood (use its NURBS surface); meshes are shifted from the NURBS/label frame (apply the shift once).
- Never truncate imaged arteries to "let growth define them" (stalls below ~60 mm); fill far tissue with XCAT.fill_supply_gaps.
- grow() falls back to root_caps = 1.588 mm (a coronary value): always set the imaged inlet radius.
- Size-dependent growth constants (connection cap, stub reach, search grid) scale with size_scale(spec) and must be numerical guards, not physiology.
- Total occlusion (area 0) -> NaN; use a tiny area fraction (1e-3).
- Body.Injection() default is 50 mL + 40 mL saline, NOT the lab protocol (0.5 mL/kg contrast at 5 mL/s, then 0.5 mL/kg of 30:70 contrast:saline); always pass it explicitly.
- Lower arterial peak at stress is real dilution (CO x1.71), not a bug.
- BasisSimulator takes <= 64 materials: DynamicPhantom.materials_at; load slabs with zrange.
- Notebook traps: no docstring before begin/let; no one-line do...end; Slider is in both PlutoUI and CairoMakie; map does not take pairs(NamedTuple); restart the worker after a package change; never pkill workers by name; CairoMakie needs XDG_CACHE_HOME.
- Docs: a new notebook needs NOTEBOOK_META + page title + committed render; sliders go in FORCE_FALLBACK_BONDS; no machine paths (/Users, /Volumes) or errors in outputs; any src/ change re-renders every page.
- Public repo: no XCAT data, no XCAT-derived files, and no wording about the lab's private perfusion-analysis code.

OPEN WORK
- Body slots and flow share for organs other than heart and lungs (brain, liver, leg).
- Two blood supplies (liver: hepatic artery + portal vein).
- Flow on very large trees (the leg: ~117 M segments) may need a hierarchical solve per patch.