PUNCH Estimates Coronary Flow Reserve Without a Wire

Analysis by the aitrendblend editorial team · AI for medical imaging and healthcare · Source paper peer reviewed in Medical Image Analysis, 2027 · Explains published research, not medical advice · 29 September 2026 · Reading time about 30 minutes

• Coronary flow reserve
• Microvascular dysfunction
• Physics informed neural networks
• Uncertainty quantification
• Coronary angiography
• Variational inference
• PyTorch

A clean angiogram is not always a clean bill of health

Coronary angiography answers one question extremely well. Is a large artery narrowed? A catheter delivers contrast dye, an X ray camera films the dye as it fills the vessels, and the cardiologist reads the shape of the tree. The paper puts the global volume at more than 10 million procedures a year, and it reports that no obstructive lesion turns up in about 70 percent of the patients evaluated for ischemic heart disease.

That number should give anyone pause. Chest pain does not require a blocked pipe. The smallest vessels, the ones that actually feed heart muscle, can fail to widen when the heart works harder. This is coronary microvascular dysfunction, or CMD, and the paper cites the EAPCI consensus document by Kunadian and colleagues (2021) for its presence in up to 70 percent of patients with angina and non obstructive arteries. Ong and colleagues (2018) are cited for the observation that it falls disproportionately on women.

The consequences are not mild. A meta analysis by Gdowski and colleagues (2020), quoted in the discussion, links CMD to nearly four times the odds of all cause death (odds ratio 3.93, 95 percent confidence interval 2.91 to 5.30) and about five times the odds of major adverse cardiac events (odds ratio 5.16, interval 2.81 to 9.47). Yet the paper notes that specialized invasive CMD assessments happen only hundreds to a few thousand times a year, concentrated in a small number of expert centers. Compare that with roughly 4 million elective coronary angiograms a year in Europe and the United States combined.

The missing measurement is coronary flow reserve, or CFR. It is the ratio of blood flow during drug induced hyperemia to flow at rest, and a low value says the microcirculation cannot open up when asked. The reference method is thermodilution. A clinician places a wire with a thermistor in the distal artery, injects room temperature saline, and reads how the temperature dips over time (Fearon and colleagues, 2003, and Barbato and colleagues, 2004). It works, but the paper is frank that it needs precise timing, depends on injection technique and temperature stability, and adds the risk of an extra instrument in the artery.

Readers of our Medical AI coverage will know the usual angiography story, where a network finds narrowings, as in the YOLO study of coronary lesion detection, or scores lesion risk, as in AngioGraphCAD. Those methods look at the shape of the artery wall. PUNCH looks at something else, how fast the dye travels.

The authors also choose their words with care. They do not call PUNCH non invasive, because the patient still undergoes angiography and still needs a vasodilator. They call it wire free. That distinction sets the scale of the claim, and it is worth keeping in mind through everything that follows.

How dye becomes a speedometer

Here is the physical idea in one sentence. Contrast dye behaves like a tracer, so the way it advances down an artery carries the speed of the blood carrying it. Indicator dilution theory has used that logic for decades (Zierler, 1962 and 2000), and thermodilution is one member of the family. The novelty is doing it from ordinary fluoroscopy frames.

The team follows the left anterior descending artery, or LAD, the vessel most often studied in CMD testing. Along the vessel centerline they stack the image intensity from every frame into a picture called a kymograph. Distance along the vessel runs one way and time runs the other. A front of dye that moves at constant speed leaves a diagonal streak, and the streak’s geometry gives the speed, roughly the distance covered divided by the time it took.

Reading that streak by eye would be fragile. Noise, cardiac motion, low frame rates and incomplete filling all blur it. So the authors write down the physics and let a network fit it. Starting from conservation of blood and of contrast in a short stretch of vessel with cross sectional area \(A(s,t)\), they get two balance equations.

Here \(I\) is contrast intensity, \(u\) is the cross sectionally averaged blood velocity, \(D\) is an effective dispersion coefficient that soaks up mixing and unresolved three dimensional effects, and \(s\) is distance along the vessel. Combining the two, then noting that the sharp gradients sit at the leading and trailing edges of the bolus, the authors arrive at a clean one dimensional advection dispersion law.

A number in the paper deserves more attention than it gets. Using a vessel length of about 0.1 m, a dispersion coefficient near \(10^{-6}\) square meters per second and a velocity near 0.04 meters per second, the effective Péclet number lands around 4000. Advection wins by a wide margin. Diffusion only matters in thin layers at the front and back of the bolus, of thickness roughly the vessel diameter divided by the square root of the Péclet number.

Read that as a design hint. In a regime this advective, the image mostly tells you where the front is at each moment, and the front is precisely what the kymograph slope measures. Dispersion is the nuisance parameter, identifiable mostly at the edges. That reading is ours and not the authors’ wording, but it explains why a one dimensional model is a reasonable bet for the proximal and mid LAD. It also says the estimate will lean hard on how sharp and how well tracked the front is.

One more detail. The paper’s Péclet estimate uses 0.04 meters per second, while its synthetic resting velocities run from 10 to 20 centimeters per second. Faster flow only raises the Péclet number, so the advection dominated picture survives either way.

The clinical quantity comes last. CFR is the ratio of the spatially averaged hyperemic velocity to the spatially averaged resting velocity along the imaged segment of length \(L\).

Because it is a ratio, any scale factor that touches both states equally, a pixel calibration error for example, should cancel in principle. The same fact has a cost. A ratio inherits the errors of both terms, and a small denominator amplifies them, a theme that also runs through our piece on why heart masks cannot lift ejection fraction prediction, where subtraction plays the same trick.

Inside the network

The idea of pushing a differential equation into a network’s loss goes back to the physics informed neural networks of Raissi, Perdikaris and Karniadakis (Journal of Computational Physics, 2019). PUNCH borrows that machinery and adds two things the authors say earlier work lacked, a measure of confidence and a way to link the rest and hyperemic studies.

Each physiological state gets its own set of three small fully connected networks with hyperbolic tangent activations. The intensity network takes position, time and a latent variable and returns dye intensity, using three hidden layers of 64 units. The velocity network takes position and the latent variable, and the dispersion network takes position, time and the latent variable. Both use two hidden layers of 32 units. Sigmoid transforms keep velocity and dispersion inside physiologically plausible ranges, and the same bounds apply to both states.

The two states, rest and hyperemia, are not independent. They share one latent variable \(\mathbf{z}\), which follows a Gaussian variational posterior learned from the data, in the spirit of the variational autoencoder of Kingma and Welling (arXiv 1312.6114). The authors describe it as an abstract per case variable that absorbs image quality, injection variability and acquisition noise. Poor data in one state therefore widens the uncertainty of the other, and the CFR uncertainty comes out coupled.

Three terms pull in different directions. The data term forces the intensity network to match the observed kymograph. The PDE term forces its derivatives, computed by automatic differentiation, to satisfy the advection dispersion law with the velocity and dispersion the other two networks propose. The KL term keeps the latent posterior from drifting away from a Gaussian prior. Each iteration draws 2000 collocation points per state for the physics residual, 2000 data points per state for the data term, and one reparameterized latent sample. Training uses Adam at an initial learning rate of \(10^{-3}\) with a cosine annealing schedule.

At inference the authors draw 100 samples of \(\mathbf{z}\), compute a CFR for each, and report the mean with a 95 percent credible interval. On a single NVIDIA A40 graphics card the whole pipeline takes about three minutes per patient.

The hyperparameters, and where they came from

The grid was modest. Prior standard deviations of 0.5, 1, 2, 5, 10 and 15 were crossed with KL weights of \(10^{-4}\), \(10^{-3}\) and \(10^{-2}\), and the full results sit in tables 6 to 22 of the supplement. The clinical evaluation used a prior standard deviation of 2, a KL weight of \(10^{-3}\) and a one dimensional latent space. That choice was made from the synthetic grid, judged by mean absolute error together with coverage and interval sharpness, and not from the patient results.

Two robustness checks stand out. Across ten random seeds on ten matched synthetic cases, the inferred CFR had a median coefficient of variation of 0.8 percent, with a mean of 1.1 percent and a worst case of 3.8 percent. Changing the training horizon among 7000, 10,000 and 13,000 epochs moved the CFR by a median of 0.6 percent between the first two and 0.1 percent between the last two. Raising the latent dimension past one or two made things worse, which the authors read as overparameterization of the latent space.

Why every patient gets a fresh model

PUNCH is trained from scratch on each case, with no parameter sharing between patients. That is unusual in deep learning and it is deliberate. There is no training set to leak into a test set, and no worry that a network trained on one population will misjudge another. The price is time, three minutes each, and the absence of any learned prior about what a typical LAD looks like.

It also changes what the word validation means. A conventional model earns trust by holding up on unseen patients. Here every patient is unseen by construction, and trust has to come from the physics being right and from checks like the ones we look at next.

The synthetic proving ground

Invasive ground truth for CMD is scarce, so the team built one. They solved the advection dispersion equation numerically, with a method of lines finite volume scheme and an explicit Runge Kutta integrator, to create 1000 synthetic angiographic sequences. Resting velocities ran from 10 to 20 centimeters per second for typical cases, and a quarter of the sequences were labeled atypical stress tests with resting velocities of 20 to 40 centimeters per second. Hyperemic flow came from scaling the resting field by a sampled CFR, drawn between 1.5 and 4.0 for typical cases and between 1.0 and 2.0 for low CFR cases. Vessel length ranged from 6 to 10 cm and simulated recordings lasted 3 to 6 seconds.

Clean simulations would flatter any method, so each kymograph was corrupted with Poisson shot noise, spatial blur, low frequency intensity drift, banding, localized dropout and temporal downsampling, over a range from mild to severe. To check the simulated CFR values were realistic, the authors compared them with the invasive CFR distribution from their 20 patients using a two sample Kolmogorov Smirnov test. The statistic was 0.19 with a p value of 0.45, so the two were statistically indistinguishable, though with 20 clinical values that test has little power to notice a difference.

On the matched synthetic set, point accuracy was excellent in the middle of the hyperparameter grid, with mean absolute errors as low as 0.035 and root mean squared errors under 0.10 for moderate prior variance and KL weight. Deming regression slopes sat between 0.96 and 0.98. At the CFR 2.5 cutoff, the posterior mean kept both sensitivity and specificity above 0.95 in most configurations, and the credible bounds held high specificity with sensitivity above 0.90 for moderate regularization.

The inverse crime, and how they answered it

There is an obvious objection. If the synthetic data come from the very equation the network enforces, of course it recovers the answer. Numerical analysts call this the inverse crime, and the authors name it themselves. Their response was to build two synthetic sets whose forward physics is richer than the inference equation. One adds a tapering vessel with a 20 to 40 percent area reduction over the segment, using the full mass conservation form. The other adds tapering plus two to four side branches that drain part of the incoming flow. PUNCH was then run unchanged.

Table 1. CFR recovery under model mismatch, from Table 1 of the paper. Lower error and a slope closer to 1 are better.

The pattern is instructive. Point estimates survive, with slopes near 0.97 and correlations of 0.98 or higher. What erodes is interval coverage, which falls from 0.672 to 0.535 and 0.574. The authors caution that the two mismatch sets are independent random draws, so the lower error with side branches does not mean that more mismatch is easier. Model error therefore hurts the honesty of the uncertainty more than it hurts the estimate. One mismatch remains untested, a hyperemic state built from a microvascular resistance model instead of a simple velocity scaling.

What the ablation really shows

Ablation studies often decorate a paper. This one carries the argument. On the first 256 synthetic cases at the clinical operating point, the team removed or replaced one component at a time.

Table 2. Ablation and uncertainty method comparison on the first 256 synthetic cases, from Table 2 of the paper. The full model coverage is computed on the 248 case subset shared by all variants.

Start with the second row. Without the PDE residual the mean absolute error explodes to 1.380 and the Deming slope collapses to 0.049. The reason is structural. The velocity and dispersion networks never touch the pixels. They receive training signal only through the physics residual, so with that term gone they are unconstrained, and the intensity network can fit the image perfectly while saying nothing about flow. It is the machine learning version of a camera that reproduces the scene and has no idea what moved.

Now the third row, which is easy to miss. A deterministic variant with no latent and no KL term is slightly more accurate than full PUNCH, with a mean absolute error of 0.021 against 0.032. Its coverage is exactly zero, because it has no interval to cover anything. The variational machinery buys uncertainty and pays a little accuracy for it. The authors say this plainly, describing the latent as a separate contribution to uncertainty, distinct from the physics that drives accuracy.

Sharing the latent across rest and hyperemia barely moves the point estimate. Independent latents give a mean absolute error of 0.034 and higher coverage of 0.746, but the intervals are wider, with a mean width of 0.095 against 0.070, so the gain comes from width and not from better calibration. The shared latent, in the authors’ argument, supplies a single source of coupled rest and hyperemia uncertainty.

The last row explains why they did not take the cheaper road. Monte Carlo dropout, in the style of Gal and Ghahramani (2016), produced a mean 95 percent interval width of 0.001 against 0.070 for the variational posterior, and coverage of 0.004. The authors trace this to structure. Per unit dropout noise averages out over the spatial domain, while the shared latent perturbs the whole velocity field coherently, which is what actually moves a spatially averaged CFR.

Twenty patients

The clinical data come from routine care at the University of Michigan Health System, analyzed retrospectively. Every patient had undergone CFR testing for a clinical reason. Angiograms were taken in matched right anterior oblique and cranial projections at 10 or 15 frames per second, once at rest and once under pharmacologic hyperemia. The invasive reference used bolus thermodilution. Three injections of 3 mL of room temperature saline were made at rest, adenosine was infused at 140 micrograms per kilogram per minute, and three more injections followed. CFR was the mean resting transit time divided by the mean hyperemic transit time.

The image pipeline has five steps. A modified segmentation network built on SAM2 with dynamic snake convolution layers (Xiong and colleagues, 2026, and Qi and colleagues, 2023) extracts the LAD centerline. LocoTrack (Cho and colleagues, 2024) follows points along that centerline from frame to frame. Multi scale Hessian filtering strengthens the dye and suppresses background structure. DICOM metadata convert pixels to millimeters. Finally the filtered frames are stacked into a kymograph. The authors verified the centerline by eye for each patient but did no manual editing, and they report that moving the start frame, end frame or sequence length within reason, as long as the window still brackets the passage of contrast, leaves the inferred CFR essentially unchanged.

Because thermodilution is itself noisy, the invasive value also got an uncertainty estimate. With three repeated measurements at rest and three during hyperemia, all 27 possible replacement combinations per state were enumerated, giving 729 bootstrap CFR values per patient.

The agreement is strong on its face. PUNCH and invasive CFR had a Spearman rank correlation of 0.89 (\(p < 10^{-6}\)), and ordinary least squares regression gave an \(R^{2}\) of 0.79. Deming regression, which allows error in both variables, produced a slope of 0.83 and an intercept of 0.11 (Linnet, 1993). A slope below one with a positive intercept means PUNCH runs a little low at the high end and sits close to the invasive value at the low end. The Bland and Altman analysis (1986) found a mean difference of minus 0.31 with a standard deviation of 0.48 and limits of agreement from minus 1.25 to 0.63. All 20 patients fell inside those limits. In relative terms the mean bias was minus 11.4 percent with a standard deviation of 20.4 percent, and 10 of 20 patients, or 50 percent, landed within 20 percent of the invasive value.

The bias has a possible partner in the reference standard. The authors point to studies, including Gallinoro and colleagues (2023), showing that bolus thermodilution tends to report higher CFR than continuous thermodilution, and suggest this may explain part of the negative offset. That is a hypothesis. This study has no continuous thermodilution measurements to test it.

The 2.5 cutoff

A clinician often wants a yes or no, so the authors also classified patients at the usual cutoff, with CFR below 2.5 counted as abnormal. Against the invasive reference, PUNCH found all 11 abnormal patients, flagged 2 normal patients as abnormal, and correctly cleared 7.

Table 3. Diagnostic performance for CMD classification at CFR below 2.5, from Table 3 of the paper. Intervals come from a 2000 sample case bootstrap.

Look closely at the two errors. Both were false positives, and both came from the negative bias. The patients with invasive CFR values of 2.8 and 3.06 were estimated at 1.93 and 2.39, just under the cutoff. Nothing here suggests PUNCH missed a sick patient, and with 11 abnormal cases that is a small sample from which to promise it never will.

The confidence intervals also deserve a raised eyebrow. Sensitivity and negative predictive value carry intervals of exactly 1.00 to 1.00, which looks like certainty. It is an artifact. A bootstrap can only resample the cases it has seen, and when none of 11 abnormal patients was missed it cannot imagine a miss. An exact binomial interval for 11 out of 11 would be far wider. The paper handles this responsibly, calling the results feasibility scale evidence and not established diagnostic performance.

A baseline that fails, and what that means

To ask whether the kymograph alone does the work, the authors built a naive image based comparator. It applies a Hough transform to find the slope of the contrast front in each kymograph and estimates CFR as the ratio of hyperemic to resting slopes, with no physics and no uncertainty.

Table 4. Image based baseline against invasive CFR, from Table 4 of the paper. Baseline metrics use only the 10 patients with a finite, positive ratio.

The Hough tracker returned a usable CFR for only half the cohort. In the other ten the discretized slope was infinite, for a near vertical front, or non physical, for a negative ratio. Where it did work, it showed no association with the invasive measurement, and the mean absolute error was 2.26 against 0.47 for PUNCH on all 20. The authors conclude that the agreement comes from the physics constrained inference and not the kymograph representation.

That is a fair reading of that comparison. It is also a comparison against a weak opponent, and the authors say why they stopped there. Frame counting methods need a standardized injection protocol that this retrospective dataset lacks, and a supervised deep learning comparator cannot be trained because no labeled angiographic CFR dataset of useful size exists. Fair enough. The consequence is that we do not yet know how PUNCH compares with the strongest angiography only method a cardiologist could build today.

Reading the uncertainty honestly

The clinical pitch of PUNCH is that each estimate arrives with a credible interval, narrow when the dye dynamics are clean and wide when the image is ambiguous. The paper proposes using width as a triage signal. A low CFR with a tight interval might justify further investigation, while a wide interval might suggest trying another projection, repeating the imaging or turning to another physiological test. That is the authors’ proposal, and the honest question is whether the intervals earn it.

On synthetic data, the answer is only partly. At the clinical operating point the raw 95 percent intervals covered the truth in 67.2 percent of matched cases, and no configuration in the hyperparameter grid reached nominal coverage. The reliability diagram in the paper shows empirical coverage below the nominal line at every level, with an expected calibration error of 0.170. The authors call these intervals under dispersed and say so in the abstract.

They then tried a repair. A post hoc split conformal step, fitted on one synthetic split and checked on a disjoint one, raised test coverage from 0.68 to 0.94, using a multiplicative factor roughly four times the nominal 1.96. That is a useful proof that the intervals can be calibrated in domain. It is not evidence that they are calibrated on patients, because the recalibration was never fitted to clinical data.

The intervals still behave sensibly in one respect. When the authors cut the effective frame rate from 15 to 7.5 frames per second, the mean absolute error rose from 0.034 to 0.086 and the root mean squared error from 0.068 to 0.158. The intervals widened roughly in proportion to the error, so coverage stayed near 0.67 and the model did not become more confident as information vanished. The sensitivity of the lower credible bound dropped from 0.95 to 0.86. The paper is careful to add that clinical acquisitions at 10 to 15 frames per second sit inside the tested range, so this is a characterization and not proof of frame rate invariance.

“PUNCH is more precise than the reference rather than disagreeing with it.”Thakur and colleagues, Medical Image Analysis, 2027

The clinical intervals tell a stranger story. PUNCH’s median 95 percent half width was about 0.04 CFR units, which is sharper than the invasive reference’s own noise. The invasive value fell inside the PUNCH interval in only 3 of 20 patients. In the other direction, the PUNCH estimate fell inside the invasive 95 percent test retest interval in 15 of 20. The Bland and Altman limits of agreement, about plus or minus 38 percent of the mean CFR, are comparable to the documented reproducibility of bolus thermodilution, roughly plus or minus 33 percent.

Put those numbers side by side. The estimates agree with the reference about as well as the reference agrees with itself. The intervals do not, because they measure something narrower. In the authors’ own words, they capture measurement and identifiability uncertainty, meaning how well the image pins down the flow, propagated through the shared latent. They do not capture the biological variability of CFR, and they do not capture doubt about whether a one dimensional model is the right model. The authors say they do not claim calibrated clinical credible intervals.

We think that framing is the right one, and it is worth internalizing. A credible interval is only as honest as the story it tells about where the noise comes from. Here the story is image noise, and the interval is best used as a quality control signal that ranks studies by reliability. The paper’s own discussion puts it plainly, calling honest AI a design principle and not a performance claim.

“the principle is partly aspirational pending the recalibration and external validation described below.”Thakur and colleagues, Medical Image Analysis, 2027

The clinical translation gap

Between a promising feasibility study and a tool a cardiologist relies on sits a long road. This paper marks the first stretch and says so. Here are the gaps we would want closed, most of them named by the authors.

The first is protocol. PUNCH removes the wire and nothing else. It still needs matched resting and hyperemic acquisitions with a vasodilator, so an angiographic procedure and a pharmacologic protocol remain. The paper states that it does not displace the need for either. Its value is therefore for centers where invasive microvascular testing is not routine, and for the patients who would otherwise never receive a functional diagnosis.

The second is anatomy. The one dimensional assumption suits the proximal and mid LAD. Severe bifurcations, vessel overlap or extreme tortuosity may break it. Side branch outflow and vessel tapering are not modeled at inference, and the synthetic mismatch tests showed they degrade interval coverage more than point accuracy. The rest and hyperemic frames are matched by projection but not by cardiac phase, so a residual mismatch is absorbed into the latent variable. ECG gated alignment is listed as future work. Other territories, the right coronary and left circumflex arteries, were not tested.

The third is the decision threshold. Two of the 20 patients sat on the wrong side of 2.5 because of the negative bias, and the authors state that a recalibration step would be advisable before any threshold based clinical use. Borderline cases near a cutoff are where a biased estimator does the most harm.

The fourth is the reference itself. Bolus thermodilution has substantial intra patient variability, with a median standard deviation of about 12 percent across repeated injections in the paper. Validating against it caps how sharply anyone can measure agreement.

Regulatory and safety notes

The paper reports no regulatory clearance and proposes no immediate clinical deployment. Software that estimates a physiological measure to support a diagnosis would normally face medical device review in most jurisdictions, and that process would ask for the multi center, multi territory validation the authors say is still needed. Readers should also note the ownership. The per patient solver belongs to AngioInsight, Inc., which funded the work, and it is not released as open source. Code that generates the synthetic dataset, including the forward solver and the mismatch generators, can be shared on reasonable request, which lets outsiders reproduce the synthetic, mismatch, ablation and uncertainty analyses. The paper also supplies pseudocode, the architecture and the operating point hyperparameters. Clinical angiograms cannot be shared for privacy reasons.

The competing interest statement is explicit. Several authors are employees of AngioInsight and hold equity in it, and others hold equity as consultants or co founders. That does not make the results wrong. It is standard practice in medical technology, and the paper discloses it. It does mean independent replication, when it comes, carries extra weight.

Where PUNCH sits among angiography derived methods

The authors place PUNCH along five axes, the imaging input, whether labeled training data are needed, run time, the physiological quantity produced, and whether uncertainty is propagated. Their Table 5 lines up representative approaches.

Table 5. Positioning among representative approaches, condensed from Table 5 of the paper.

Two things stand out. PUNCH is the only entry that propagates uncertainty, and it is the slowest of the label free angiography methods, at minutes where the others take seconds. Computational fluid dynamics methods need high resolution CT and hours of computation per patient (Taylor and colleagues, 2013), and supervised approaches need labeled data that for CMD phenotypes essentially do not exist. PUNCH sidesteps both by learning from a single patient’s own images. A ratio of flows also asks a different question than a pressure ratio like FFR, since it speaks to the microcirculation and not only to a narrowing.

Where the idea could travel

The recipe is portable. Take a governing equation, fit it per case rather than across a population, and add a shared latent variable so every estimate arrives with a statement of doubt. The authors say this offers an alternative to label hungry supervised learning whenever ground truth is scarce, and they name perfusion imaging, tracer kinetics and dynamic contrast studies in other vascular beds as candidates. They are equally clear that whether it extends is an open question they do not resolve.

Our sense is that two neighbors on this site frame the trade off well. A neural operator for free boundary problems learns a solver that generalizes across problem instances, the opposite bet from PUNCH’s fresh network per patient. And the physics guided synthetic ultrasound study shows the same appetite for simulators that generate labeled data where clinics cannot. Operator learning amortizes cost and needs a training distribution. Per case fitting pays each time and needs none. Which wins depends on how much you trust the population you trained on.

Limitations, sample size and bias

Set the achievements aside for a moment and read the constraints in the paper’s own numbers.

• Sample size. The clinical cohort is 20 patients, with 11 abnormal and 9 normal by the invasive reference. Specificity has a confidence interval from 0.50 to 1.00, and the sensitivity interval of 1.00 to 1.00 is a bootstrap artifact.
• Dataset bias. One center, one artery, retrospective data, and patients who were referred for CFR testing for clinical reasons. Such a cohort is not a screening population, and its spread of CFR values may differ from the patients most in need.
• Circular evidence. The synthetic data come from the same family of equations that PUNCH enforces at inference. The mismatch experiments reduce that concern without removing it, and the hyperparameters were chosen on the synthetic grid.
• Calibration. Raw intervals cover 67.2 percent where 95 percent is nominal, the recalibration was validated on synthetic data only, and the intervals are sharper than the invasive reference’s own variability.
• Systematic bias. A mean difference of minus 0.31 pushed two patients across the diagnostic cutoff.
• Generalization. Tested at 10 to 15 frames per second and in matched right anterior oblique and cranial projections, with clear degradation at 7.5 frames per second. Other frame rates, other views and other coronary territories are untested.
• Baselines. The comparison covers a naive Hough baseline only, with no frame counting or supervised comparator.
• Transparency. The clinical data cannot be shared and the per patient solver is proprietary, so outside groups cannot rerun the central result today.

None of these is hidden. Almost every item above is stated in the paper’s own limitations section, and the authors go as far as calling their principle of honest uncertainty partly aspirational. That candor is the strongest reason to take the work seriously as a first step.

A PyTorch reconstruction of PUNCH

The authors did not release their per patient solver, so there is no official repository to point to. What follows is our own reconstruction from the paper’s equations, network sizes and training protocol, written so you can see every moving part. Where the paper leaves a detail open, such as initialization, output heads, the log scale bound on dispersion and the resolution of the toy simulator, the comments say DESIGN CHOICE. This is educational code and not clinical software, and it has never touched patient data.

The file contains the building blocks, the two branch model with a shared variational latent, the three part loss and the per case training loop, Monte Carlo inference of CFR with a credible interval, cohort evaluation helpers, a split conformal recalibration factor, a toy advection dispersion simulator, and a runnable smoke test on dummy data. It needs Python 3, PyTorch 2 and NumPy.

"""
PUNCH reference sketch. Physics informed, uncertainty aware estimation of coronary flow reserve (CFR).

Reconstructed from the description in
    Thakur et al. PUNCH, Physics informed uncertainty aware network for coronary hemodynamics.
    Medical Image Analysis 115 (2027) 104273.  https://doi.org/10.1016/j.media.2026.104273

IMPORTANT. The authors' per patient solver is proprietary and was not released. This file is an
independent educational reconstruction from the paper's equations, network sizes and training
protocol. Details the paper does not state (initialisation, output heads, the log scale mapping
for dispersion, the grid resolution of the toy simulator) are our own choices and are marked
with "DESIGN CHOICE". It is not clinical software and it has not been validated on patient data.

Layout
    1. Building blocks (MLP, bounded outputs, one branch per physiological state)
    2. The PUNCH model with a shared variational latent
    3. Losses (data fidelity, advection dispersion residual, KL) and the per case training loop
    4. Monte Carlo inference of CFR with a 95 percent credible interval
    5. Evaluation helpers and a split conformal recalibration factor
    6. A toy advection dispersion simulator and a runnable smoke test on dummy data
"""
import math
import numpy as np
import torch
import torch.nn as nn


# ---------------------------------------------------------------------------
# 1. Building blocks
# ---------------------------------------------------------------------------
class MLP(nn.Module):
    """Fully connected network with hyperbolic tangent activations, as in the paper."""

    def __init__(self, in_dim, hidden, out_dim=1):
        super().__init__()
        dims = [in_dim] + list(hidden)
        layers = []
        for a, b in zip(dims[:-1], dims[1:]):
            layers += [nn.Linear(a, b), nn.Tanh()]
        layers.append(nn.Linear(dims[-1], out_dim))
        self.net = nn.Sequential(*layers)

    def forward(self, x):
        return self.net(x)


def bounded_linear(raw, lo, hi):
    """Sigmoid parameterisation that keeps an output inside [lo, hi]."""
    return lo + (hi - lo) * torch.sigmoid(raw)


def bounded_log10(raw, lo, hi):
    """DESIGN CHOICE. Sigmoid on a log10 scale so that D can span an order of magnitude range."""
    return 10.0 ** (math.log10(lo) + (math.log10(hi) - math.log10(lo)) * torch.sigmoid(raw))


class Branch(nn.Module):
    """One physiological state (rest or hyperaemia) with three subnetworks.

    intensity network I(s, t; z)   three hidden layers of 64 units
    velocity network  u(s; z)      two hidden layers of 32 units
    dispersion network D(s, t; z)  two hidden layers of 32 units
    """

    def __init__(self, z_dim=1, u_bounds=(0.2, 20.0), d_bounds=(1e-5, 1e-1)):
        super().__init__()
        self.net_i = MLP(2 + z_dim, (64, 64, 64))
        self.net_u = MLP(1 + z_dim, (32, 32))
        self.net_d = MLP(2 + z_dim, (32, 32))
        self.u_bounds, self.d_bounds = u_bounds, d_bounds

    def intensity(self, s, t, z):
        return self.net_i(torch.cat([s, t, z.expand(s.shape[0], -1)], dim=1))

    def velocity(self, s, z):
        raw = self.net_u(torch.cat([s, z.expand(s.shape[0], -1)], dim=1))
        return bounded_linear(raw, *self.u_bounds)

    def dispersion(self, s, t, z):
        raw = self.net_d(torch.cat([s, t, z.expand(s.shape[0], -1)], dim=1))
        return bounded_log10(raw, *self.d_bounds)


# ---------------------------------------------------------------------------
# 2. The PUNCH model. Two branches that share one variational latent z
# ---------------------------------------------------------------------------
class PUNCH(nn.Module):
    def __init__(self, z_dim=1, prior_mean=45.0, prior_std=2.0, u_bounds=(0.2, 20.0), d_bounds=(1e-5, 1e-1)):
        super().__init__()
        self.z_dim = z_dim
        self.rest = Branch(z_dim, u_bounds, d_bounds)
        self.hyper = Branch(z_dim, u_bounds, d_bounds)
        # Variational posterior q(z) = N(mu, diag(sigma^2)). Free parameters, one set per case.
        # DESIGN CHOICE. Start the posterior at the prior.
        self.mu_z = nn.Parameter(torch.full((z_dim,), float(prior_mean)))
        self.logvar_z = nn.Parameter(torch.full((z_dim,), math.log(prior_std ** 2)))
        self.prior_mean, self.prior_std = prior_mean, prior_std

    def sample_z(self, n=1):
        """Reparameterised sample z = mu + sigma * eps, shape (n, z_dim)."""
        eps = torch.randn(n, self.z_dim, device=self.mu_z.device)
        return self.mu_z + torch.exp(0.5 * self.logvar_z) * eps

    def kl_to_prior(self):
        """Closed form KL( N(mu, sigma^2) || N(mu_p, sigma_p^2) ), summed over latent dimensions."""
        var, var_p = torch.exp(self.logvar_z), self.prior_std ** 2
        kl = 0.5 * (torch.log(var_p / var) + (var + (self.mu_z - self.prior_mean) ** 2) / var_p - 1.0)
        return kl.sum()


# ---------------------------------------------------------------------------
# 3. Losses and per case training
# ---------------------------------------------------------------------------
def pde_residual(branch, s, t, z):
    """Advection dispersion residual  I_t + u I_s - D I_ss  from Eq. 3, by automatic differentiation."""
    s = s.clone().requires_grad_(True)
    t = t.clone().requires_grad_(True)
    i = branch.intensity(s, t, z)
    ones = torch.ones_like(i)
    i_s = torch.autograd.grad(i, s, ones, create_graph=True)[0]
    i_t = torch.autograd.grad(i, t, ones, create_graph=True)[0]
    i_ss = torch.autograd.grad(i_s, s, torch.ones_like(i_s), create_graph=True)[0]
    u = branch.velocity(s, z)
    d = branch.dispersion(s, t, z)
    return i_t + u * i_s - d * i_ss


def kymograph_points(kymo, n_points, generator=None):
    """Sample (s, t, intensity) triplets uniformly from a kymograph of shape (n_t, n_s).

    Coordinates are normalised to the unit square, s by arc length and t by recording duration.
    Intensities are normalised to the unit interval, as the paper describes.
    """
    kymo = np.asarray(kymo, dtype=np.float32)
    kymo = (kymo - kymo.min()) / (kymo.max() - kymo.min() + 1e-8)
    n_t, n_s = kymo.shape
    it = torch.randint(0, n_t, (n_points,), generator=generator)
    js = torch.randint(0, n_s, (n_points,), generator=generator)
    s = (js.float() / (n_s - 1)).unsqueeze(1)
    t = (it.float() / (n_t - 1)).unsqueeze(1)
    val = torch.from_numpy(kymo)[it, js].unsqueeze(1)
    return s, t, val


def punch_loss(model, rest_kymo, hyper_kymo, n_f=2000, n_d=2000, kl_weight=1e-3):
    """L = L_data + L_PDE + lambda_KL * KL(q(z) || p(z)), Eq. 4, summed over both states."""
    z = model.sample_z(1)  # one reparameterised latent sample per iteration
    parts = {"data": 0.0, "pde": 0.0}
    total = 0.0
    for branch, kymo in ((model.rest, rest_kymo), (model.hyper, hyper_kymo)):
        s_d, t_d, i_d = kymograph_points(kymo, n_d)
        l_data = torch.mean((branch.intensity(s_d, t_d, z) - i_d) ** 2)
        s_f, t_f = torch.rand(n_f, 1), torch.rand(n_f, 1)  # collocation points, uniform on (s, t)
        l_pde = torch.mean(pde_residual(branch, s_f, t_f, z) ** 2)
        total = total + l_data + l_pde
        parts["data"] += float(l_data.detach())
        parts["pde"] += float(l_pde.detach())
    kl = model.kl_to_prior()
    total = total + kl_weight * kl
    parts["kl"] = float(kl.detach())
    return total, parts


def fit_case(rest_kymo, hyper_kymo, epochs=10000, lr=1e-3, n_f=2000, n_d=2000, kl_weight=1e-3,
             prior_mean=45.0, prior_std=2.0, z_dim=1, seed=0, log_every=0):
    """Train PUNCH independently for ONE patient. No parameters are shared across patients.

    DESIGN CHOICE. The default of 10000 epochs is the middle of the 7000, 10000 and 13000 horizons
    the paper compares. The paper does not tie the default to a single number."""
    torch.manual_seed(seed)
    model = PUNCH(z_dim=z_dim, prior_mean=prior_mean, prior_std=prior_std)
    opt = torch.optim.Adam(model.parameters(), lr=lr)
    sched = torch.optim.lr_scheduler.CosineAnnealingLR(opt, T_max=epochs)  # cosine annealing
    history = []
    for ep in range(epochs):
        opt.zero_grad()
        loss, parts = punch_loss(model, rest_kymo, hyper_kymo, n_f, n_d, kl_weight)
        loss.backward()
        opt.step()
        sched.step()
        history.append(float(loss.detach()))
        if log_every and (ep % log_every == 0 or ep == epochs - 1):
            print(f"  epoch {ep:6d}  loss {float(loss.detach()):.5f}  data {parts['data']:.5f}  "
                  f"pde {parts['pde']:.5f}  kl {parts['kl']:.3f}")
    return model, history


# ---------------------------------------------------------------------------
# 4. Monte Carlo inference of CFR and its credible interval
# ---------------------------------------------------------------------------
@torch.no_grad()
def infer_cfr(model, duration_rest_s, duration_hyper_s, n_samples=100, n_grid=200):
    """Draw N posterior samples of z, compute CFR_i = mean(u_hyper) / mean(u_rest), report the
    posterior mean and the percentile 95 percent credible interval.

    Velocities are learned in normalised units (arc length 1, duration 1). Converting to physical
    units multiplies each state by (length / duration). The length cancels in the ratio, so
    CFR = (u_h_norm / u_r_norm) * (T_rest / T_hyper).
    """
    s = torch.linspace(0, 1, n_grid).unsqueeze(1)
    cfr = []
    for _ in range(n_samples):
        z = model.sample_z(1)
        u_r = model.rest.velocity(s, z).mean()   # spatial mean over the vessel, Eq. 5 numerator/denominator
        u_h = model.hyper.velocity(s, z).mean()
        cfr.append(float(u_h / u_r) * (duration_rest_s / duration_hyper_s))
    cfr = np.array(cfr)
    lo, hi = np.percentile(cfr, [2.5, 97.5])
    return {"mean": float(cfr.mean()), "lo": float(lo), "hi": float(hi), "samples": cfr}


# ---------------------------------------------------------------------------
# 5. Evaluation helpers and split conformal recalibration
# ---------------------------------------------------------------------------
def _rank(x):
    return np.argsort(np.argsort(x)).astype(float)


def evaluate_cohort(mean, lo, hi, truth, cutoff=2.5):
    """Point accuracy, agreement, interval calibration and the CFR < cutoff classification."""
    mean, lo, hi, truth = map(np.asarray, (mean, lo, hi, truth))
    err = mean - truth
    diff_sd = err.std(ddof=1) if len(err) > 1 else float("nan")
    pred_pos, true_pos = mean < cutoff, truth < cutoff
    tp, fp = int((pred_pos & true_pos).sum()), int((pred_pos & ~true_pos).sum())
    fn, tn = int((~pred_pos & true_pos).sum()), int((~pred_pos & ~true_pos).sum())
    return {
        "MAE": float(np.abs(err).mean()),
        "RMSE": float(np.sqrt((err ** 2).mean())),
        "bias": float(err.mean()),
        "pearson_r": float(np.corrcoef(mean, truth)[0, 1]) if len(mean) > 2 else float("nan"),
        "spearman_rho": float(np.corrcoef(_rank(mean), _rank(truth))[0, 1]) if len(mean) > 2 else float("nan"),
        "bland_altman_limits": (float(err.mean() - 1.96 * diff_sd), float(err.mean() + 1.96 * diff_sd)),
        "coverage_95": float(((truth >= lo) & (truth <= hi)).mean()),
        "mean_interval_width": float((hi - lo).mean()),
        "confusion_TP_FP_FN_TN": (tp, fp, fn, tn),
        "sensitivity": tp / max(tp + fn, 1),
        "specificity": tn / max(tn + fp, 1),
    }


def fit_conformal_scale(mean, lo, hi, truth, alpha=0.05):
    """Split conformal factor k for the half width, fitted on a held out calibration set.

    Score_i = |truth_i - mean_i| / halfwidth_i. The corrected interval is mean +/- k * halfwidth.
    This follows the recalibration idea described in the paper's Limitations, not its exact code.
    """
    mean, lo, hi, truth = map(np.asarray, (mean, lo, hi, truth))
    half = np.maximum((hi - lo) / 2.0, 1e-9)
    scores = np.abs(truth - mean) / half
    n = len(scores)
    q = min(1.0, math.ceil((n + 1) * (1 - alpha)) / n)
    return float(np.quantile(scores, q))


def apply_conformal_scale(mean, lo, hi, k):
    mean, lo, hi = map(np.asarray, (mean, lo, hi))
    half = (hi - lo) / 2.0
    return mean - k * half, mean + k * half


# ---------------------------------------------------------------------------
# 6. Toy simulator and smoke test
# ---------------------------------------------------------------------------
def simulate_kymograph(u_mean, length_m=0.08, duration_s=5.0, n_s=64, n_frames=60, d_coef=1e-6,
                       pulse_on_s=0.3, pulse_len_s=1.6, noise=0.05, seed=0):
    """Toy 1D advection dispersion solver. First order upwind advection, central dispersion,
    explicit stepping. Returns a noisy kymograph of shape (n_frames, n_s) and the clean one.

    DESIGN CHOICE. This is far simpler than the method of lines Runge Kutta solver in the paper.
    Numerical diffusion from the upwind scheme dominates the physical D here, which is acceptable
    for a smoke test and would not be acceptable for generating a benchmark.
    """
    rng = np.random.default_rng(seed)
    ds = length_m / n_s
    u = u_mean * (1.0 + 0.15 * np.sin(2 * np.pi * np.linspace(0, 1, n_s)))   # gentle velocity variation
    dt = 0.4 * ds / u.max()
    n_steps = int(duration_s / dt)
    frame_every = max(n_steps // n_frames, 1)
    c = np.zeros(n_s)
    frames = []
    for k in range(n_steps):
        t = k * dt
        inlet = 1.0 if pulse_on_s <= t <= pulse_on_s + pulse_len_s else 0.0
        c_up = np.concatenate([[inlet], c[:-1]])
        adv = -u * (c - c_up) / ds
        disp = d_coef * (np.concatenate([[inlet], c[:-1]]) - 2 * c + np.concatenate([c[1:], [c[-1]]])) / ds ** 2
        c = c + dt * (adv + disp)
        if k % frame_every == 0 and len(frames) < n_frames:
            frames.append(c.copy())
    clean = np.array(frames)
    noisy = np.clip(clean + rng.normal(0, noise, clean.shape), 0, None)   # additive sensor noise
    return noisy, clean


def smoke_test(epochs=400, n_points=600):
    """Runs end to end on dummy data. It checks that every component trains and connects.
    It does NOT check clinical accuracy, and 400 epochs is far short of the thousands of epochs
    used in the paper's convergence checks (7000, 10000 and 13000)."""
    torch.manual_seed(0)
    np.random.seed(0)
    duration = 5.0
    true_cfr = 2.5
    rest, _ = simulate_kymograph(u_mean=0.04, duration_s=duration, seed=1)
    hyper, _ = simulate_kymograph(u_mean=0.04 * true_cfr, duration_s=duration, seed=2)
    print(f"kymograph shapes  rest {rest.shape}  hyper {hyper.shape}  true CFR {true_cfr}")

    model, hist = fit_case(rest, hyper, epochs=epochs, n_f=n_points, n_d=n_points, log_every=100)
    assert np.isfinite(hist).all(), "loss went non finite"
    assert hist[-1] < hist[0], "loss did not decrease"

    out = infer_cfr(model, duration, duration, n_samples=100)
    assert out["lo"] <= out["mean"] <= out["hi"], "interval does not bracket the mean"
    print(f"CFR estimate {out['mean']:.2f}  95% interval [{out['lo']:.2f}, {out['hi']:.2f}]  "
          f"(untrained accuracy, {epochs} epochs)")

    # Cohort level helpers on dummy numbers.
    rng = np.random.default_rng(0)
    truth = rng.uniform(1.0, 4.0, 40)
    mean = truth + rng.normal(-0.1, 0.3, 40)
    half = np.full(40, 0.15)                       # deliberately over confident intervals
    lo, hi = mean - half, mean + half
    metrics = evaluate_cohort(mean, lo, hi, truth)
    k = fit_conformal_scale(mean[:20], lo[:20], hi[:20], truth[:20])
    lo2, hi2 = apply_conformal_scale(mean[20:], lo[20:], hi[20:], k)
    cover_after = float(((truth[20:] >= lo2) & (truth[20:] <= hi2)).mean())
    print(f"dummy cohort  MAE {metrics['MAE']:.3f}  raw coverage {metrics['coverage_95']:.2f}  "
          f"conformal factor {k:.2f}  held out coverage after {cover_after:.2f}")
    assert k > 1.0, "over confident intervals should be widened"
    print("smoke test passed")


if __name__ == "__main__":
    smoke_test()

The smoke test is a plumbing check. It confirms that the loss is finite and falls, that the interval brackets the mean, that the cohort metrics compute and that conformal scaling widens deliberately overconfident intervals. The dummy cohort in the last step is random numbers, so its coverage figures describe the helper functions and nothing else.

The CFR estimate of 13.05 against a true value of 2.5 is wrong, and it is meant to be. Four hundred epochs on a laptop is far short of what the paper’s convergence checks used. In a separate run of 3000 epochs on the same toy case with 1000 collocation points and 1000 data points per state, our reconstruction returned 2.52 with an interval of about 2.51 to 2.52. One toy case shows that the components connect and that the physics term can recover a velocity ratio from noisy frames. It says nothing about patients, and the very narrow interval is a small reminder of the calibration problem the paper documents.

Where this leaves us

The core achievement of PUNCH is easy to state. From paired resting and hyperemic angiograms, with no wire and no training set, a small network fitted to one patient at a time produced CFR estimates that tracked invasive thermodilution at a rank correlation of 0.89 across 20 patients, in about three minutes each. It did so while attaching an interval to every estimate. The ablation adds something more valuable than the headline, a clean account of which part does what. The transport equation makes the estimate accurate, and the variational latent makes doubt expressible.

The conceptual shift matters as much as the number. Most medical imaging AI learns a mapping from pictures to labels and hopes the population it saw resembles the one it will meet. PUNCH poses a physical inverse problem for each patient, asks which velocities would make dye behave as it did, and reads a clinical quantity off the answer. Uncertainty stops being an afterthought bolted onto a classifier and becomes a designed output. The authors keep honest AI as a principle and refuse to call it an achievement, and that restraint reads as maturity.

The idea should travel, with conditions. Anywhere a trustworthy governing equation links an image to a hidden quantity and labels are scarce, the same recipe applies, with perfusion imaging and tracer kinetics as the natural next candidates. The per patient philosophy has a cousin in the cardiac digital twins we covered earlier, where a model is tuned to one heart instead of an average one. The condition is that the equation must hold in the domain, and the mismatch experiments here show that when it does not, the point estimate degrades gracefully while the uncertainty degrades badly.

The remaining limitations are not small, and the paper does not treat them as small. Twenty patients from a single center cannot establish diagnostic performance. The raw intervals are too narrow, the recalibration was only tested on synthetic data, and a bias of minus 0.31 pushed two patients over the 2.5 cutoff. The strongest comparator, a frame counting or supervised method, is missing, and outside groups cannot rerun the clinical result while the solver stays proprietary. Those facts belong next to the correlation coefficient whenever the work is discussed.

The road ahead is laid out clearly. The authors call for multi center, multi territory validation, ECG gated alignment of the rest and hyperemic runs, explicit modeling of tapering and side branches, and a microvascular resistance based construction of the hyperemic state. They also propose, without applying it yet, a composite model selection score that balances mean absolute error, coverage at 95 percent, mean interval width as a sharpness penalty, and threshold based decision performance. Fitting the recalibration to real clinical data, and testing it on new patients, would be the single most persuasive next experiment.

A wire measures flow by going into the artery. PUNCH asks the picture to say what the wire would have said, and the most encouraging thing about it is how carefully it admits when it is unsure.