This started as a final project for CS066 and did not stop when the semester did. The first version of this post was about learning the method: reimplementing Raissi, reading DeepXDE's source, building a Stiff-PINN curriculum for a combustion ODE, and running sweeps on a GPU that was dramatically overqualified for the job. That's still the first section. The rest is what happened when I pointed all of it at one real model for a summer, with Nate Kim and our advisor, and found out where it stops working.
spring
A physics-informed neural network is trained to satisfy a differential equation, not just to fit points. You move the equation to one side so a correct solution makes it zero, evaluate that on the network's output at scattered points, and add whatever's left to the loss. Autodiff gives you the derivatives of the network with respect to its inputs for free. No mesh, no discretization, and the physics does the job data usually does, which is why you can fit with very few noisy observations.
The lessons that survived: the weight on the physics term relative to the data term changes everything, adaptive weighting by gradient magnitude beat fixed weights on everything I tried, Fourier feature embeddings fix the case where a plain MLP learns the smooth part and stalls on the oscillatory part, and L-BFGS after Adam converges better than either alone. Also that the GPU is not for the network, which is tiny. It's for running thirty experiments at once.
the model
Seven ordinary differential equations for a colorectal cancer signaling network. The states are β-catenin, APC, HOXA5, HOXA13, MYC, retinoic acid, and CYP26A1, and they form three loops: a pro-stemness loop where β-catenin drives MYC drives HOXA13 drives β-catenin, a differentiation loop where retinoic acid drives HOXA5 drives APC which degrades β-catenin, and a homeostasis loop where RA induces the enzyme that clears it. The full model our advisor had was 27 species in MATLAB. The first thing I did was port it to Python, solve_ivp with a stiff solver, and match the MATLAB trajectories to under 1% error, because I didn't want to be debugging a translation for the rest of the summer. The seven-state version is a reduction of that, down to the nodes the biology question actually turns on, nondimensionalized so one time unit is roughly an hour.
Written out, with H a Hill function and d_P the APC degradation multiplier:
Every equation is production minus removal. The ε's are timescales, the η's are how strongly one thing drives another, the κ's are where each Hill curve turns over. RA gets a background input, a 24-hour circadian term, and a smooth ATRA treatment pulse; u_W and u_M are the WNT and MYC inputs, zero unless an experiment turns them on. Horizon is 150 time units. The reference trajectories come from Radau at tight tolerances, and I want to be clear that solving this system is not the hard part. A stiff solver does it in milliseconds.
Two knobs define disease. W is the WNT drive coming into the system and θP is how much functional APC you have left. Normal is (0.8, 1.0). Early adenoma is (1.0, 0.75). Advanced adenoma is (1.5, 0.5). Severe APC loss is (2.0, 0.25). Everything downstream, including a stemness score built from the states, follows from those two numbers plus 34 rate constants, Hill thresholds, and time scales. That's the 36 unknowns.
The forward problem is: given all 36, simulate. The inverse problem is: given noisy time series from an experiment, recover the 36. The inverse problem is the one that matters, because the biology question is "how much APC has this tumor lost," and that's θP.
forward
The forward PINN is a four-layer, width-256 GELU network on a Fourier embedding of time, 16 sine and 16 cosine features with frequencies drawn once and frozen. That embedding exists because the plain MLP learned the trend and missed the circadian ripple entirely: its RA output spanned 0.16 to 0.44 where the truth spanned 0.16 to 0.79. With 40 irregular observation times per regime, no noise, 50,000 collocation points, Adam then L-BFGS in float64, mean relative L2 error across all states and regimes is 2.41%. A supervised baseline with 100 evenly spaced labels and no physics gets 1.06%. So the physics term is doing something with less than half the data, but I never ran the matched 40-point physics-off control, so I can't tell you exactly how much.
One debugging story from June that changed how I look at plots. The Fourier checkpoint had converged, data MSE 1.45e-5, and I spent a day thinking it hadn't, because the PNGs on disk were from the previous plain-MLP run and I never regenerated them. The config, the checkpoint, and the figure have to come from the same run or none of them mean anything. I now write all three in the same directory with the same timestamp.
inverse
The first thing I learned this summer is that the picture above means much less than it looks like it means. The network reproduces the β-catenin curve to within the noise, and the parameters that produced that curve are off by 20% in θP and over 100% in a few of the rate constants. Lots of different parameter vectors produce the same trajectory. That isn't a PINN problem, it's the problem, and every method I tried this summer was a different way of asking how bad it is.
Early on we ran a two-parameter version, just W and θP with everything else known. W came back fine in every regime. θP came back within 1% in the normal regime and 142% off in severe APC loss. Part of the reason is in the equation: APC degradation is 1 + δ(1 − θP), so at θP = 1 the δ term vanishes and θP and δ can't be told apart at all when both are unknown. In the two-parameter run δ was fixed, so that's not the explanation there; there it's the plainer thing, that at high WNT drive the system saturates and large changes in θP barely move the curves. Either way the parameter the model is most sensitive to overall turns out to be the least identifiable where it matters. That θP wall shows up in every method below.
The full 36-parameter inverse problem, with an autodiff residual, got 4 of 36 within 10% in the normal regime and 8, 7, 7 in the others. That's the baseline.
the integral residual
I spent a while assuming that was a data problem. It wasn't. Denser, cleaner observations made W worse, not better, going from 0.80 down to 0.50. Gradient matching on frozen networks gave 0.63. Something in the loss itself was biased.
The residual in a standard PINN uses dz/dt computed by autodiff through the network. The network is a smooth function that passes close to the data, but its derivative is nobody's target and it's noisy in a way that doesn't average out. When you ask the optimizer to make dz/dt − f(z, θ) zero and dz/dt is wrong by a small systematic amount, the cheapest fix is to bend θ, not the network. You get a beautiful trajectory and a biased parameter.
The fix is to not differentiate. Between adjacent collocation points the residual becomes the trapezoid rule:
It only uses values the network outputs, not their derivatives, and it's the same integral form as the ODE solver that generated the reference data in the first place. Same networks, same optimizer, same ten conditions, and counted as parameters within 10% across the four regimes, the autodiff version had 37 and the integral version has 50. Not quite double. It's also a bundle, not a clean ablation: the same run switched to a fixed sorted collocation grid and per-state relative weighting, and the control that would have separated those out crashed. That's the one methodological result I'd defend from the summer, with that asterisk on it. A classical least-squares ODE fit with a stiff solver is still a stronger point estimator on the same data, 82 across the four regimes, and the PINN's real advantage isn't a point estimate. It's that the same machinery extends to the Bayesian and hybrid versions below, which the classical fit doesn't.
One more thing about that 50: it's the start picked by lowest physics loss out of two, and in three of four regimes that rule picked the start with worse parameter recovery. Picking by ground truth would give 55. You don't get to pick by ground truth.
The other thing that helped was experimental design, and it helped less than I expected. We started with six conditions that differ only in the retinoic acid protocol: control, no ATRA, early ATRA, late ATRA, low dose, strong circadian. Every one of them enters the WNT/HOX/MYC sub-network through a single door, the RA to HOXA5 term, so the parameters behind that door don't move. Adding WNT pulse, WNT inhibition, MYC pulse and a combo condition helped the classical fit a lot, 61 to 82, and the PINN a little, 26 to 37. Somewhere in the middle I also tried shrinking the networks from 208k weights to 38k with weight decay and alternating network and parameter updates, on the theory that the networks were too flexible. 25. Wasn't that.
identifiability
Fisher information. Finite-difference sensitivities of every state to every parameter at 200 time points across all ten conditions, 14,000 observations, FIM = SᵀS/σ² in log-parameter coordinates. The condition numbers are 10⁷ to 10¹⁹. Each regime has exactly one null direction, and it's the θP/δ direction: pure δ at normal, mixing into θP as WNT drive rises. The number of parameter pairs correlated above 0.99 goes 0, 1, 6, 10 across the regimes. Fix 28 parameters at truth and the remaining eight have no null direction at all. The first run of this script said all 36 were identifiable because I read the diagonal of a pseudoinverse; I rewrote it the same afternoon to look at null-space participation instead.
Profile likelihoods. Fix one parameter at a grid of values, refit the other 35, see whether the likelihood has a real minimum. Slow and the most honest of the four. Only the advanced adenoma regime finished. Twenty come back identifiable, sixteen weak, none flat. εM's interval is [0.596, 0.614] around a truth of 0.60. θP's is [0.16, 0.79] around 0.50. Same wall. I'd trust the shapes more than the interval endpoints; some of them sit right on the optimum, which usually means the nuisance refit didn't fully converge.
Sobol sensitivity. Global variance decomposition with SALib: sample the parameter space, simulate, and ask what fraction of the variance in each output each parameter explains, alone and with interactions. The first run said θP dominated the stemness score, first-order index 0.75, and 13 of the 14 parameters I varied tripped the nonlinearity flag. Then I noticed the ±30% box put θP as high as 1.3, where 1 + 3.5(1 − θP) goes negative and APC grows instead of decaying. Those samples were a different dynamical system. Clip the box at 1 and θP's first-order index on stemness drops to 0.15, behind η13 at 0.23, the nonlinearity flags go to zero, and the first-order indices sum to 0.95. θP still has the biggest local elasticity, about −3.6. It's just not the whole story anymore. Sensitivity and identifiability are different questions either way, and this model is the cleanest example I've seen of a parameter scoring high on one and low on the other, because the response to θP saturates.
Those three agree with each other, which is the point of running three. The fourth is the one that puts error bars on it.
Bayesian PINN. Freeze the networks from the integral-residual run, treat the 36 parameters as unknowns with log-normal priors, and run Hamiltonian Monte Carlo with the trapezoidal residual as the likelihood. The first run reported 36 of 36 identifiable with an effective sample size of three. I'd used the raw residual count, 209,930 terms, as if they were independent observations, which makes the likelihood absurdly sharp and the chain stuck. Refining the collocation grid would have made the posterior sharper without adding a single measurement. The fix is a fixed effective count, 40 per state per condition, 2,800, scaled so reduced chi-squared is one at the MAP. After that the chains mix on some parameters and not others, median ESS 96, 1265, 12, 15 across the regimes, and no run passes the gate I set for it (median ESS above 200 and credible intervals covering at least 33 of 36 truths), even at 3,000 draws and 48 leapfrog steps, which is three and a half hours per regime.
The picture is worth more than the gate. W is tight in three regimes and biased low in the fourth. θP either reverts to the prior or converges confidently to 1.0 when the truth is 0.25, a z-score in the thousands. That's the same wall as the two-parameter experiment in June, now with error bars. And it's one chain per regime, so no R-hat; I know.
The forward Bayesian PINN, parameters known and the network weights sampled, worked. Predictive coverage 0.94 to 0.95 in every regime, with the caveat that you have to report the predictive band, not the reconstruction band. The reconstruction band alone under-covers at 0.56 to 0.73 because it only carries epistemic uncertainty and ignores the observation noise. I made that mistake first and the coverage number looked like a bug for a day. The other caveat is that the 40 points it's scored on are the same 40 it was fit on. I called that held-out in an earlier draft. It isn't.
hybrids
The last thing we tried is replacing one Hill term with a small network and asking whether it learns the right shape. Universal differential equations, basically. The network is two layers of five tanh units, 46 weights, wrapped as softplus(g(x)) · (1 − e^(−x/x₀)) so it's non-negative and zero at zero, trained jointly with the remaining 34 mechanistic parameters. Two scores per run: how many of the 34 come back within 10%, and the NRMSE of the learned curve against the true Hill function over the range the data actually visited.
Replacing the WNT-to-MYC term, η b²/(κ² + b²), gives a learned curve within 6 to 10% of the true Hill function across the observed range. It also costs you the rest of the MYC equation. The basal MYC production constant comes back 32% wrong in the normal regime, against 0.4% for the mechanistic control, because the learned function absorbs it: the network can produce a constant offset at low β-catenin as easily as it can produce a curve, and nothing in the data distinguishes the two.
The zero-at-zero anchor is supposed to forbid that. It forbids a global shift. It doesn't forbid a shift over the range you observed with a bend back to zero in the range you didn't, and the data never gets close to zero. Slide the observed range around below and watch the two curves become indistinguishable.
Both functions are zero at zero.
The measurements decide if that matters.
Away from the anchor, a changed basal rate can cancel almost the entire offset. The largest remaining production mismatch in the shaded range is 0.000.
Replacing the RA-to-HOXA5 term learns the right shape 8 to 24% off depending on regime, and the anchor never binds because the gate is already at 0.7 at the lowest RA we see. The unconstrained version fits better. Recovery on the other 34 was actually a hair better than the control there, 53 to 51, mostly because the control was bad at a5 to begin with.
The one hybrid that worked cleanly is the APC one, and it worked because of how it was trained. The unknown is the degradation multiplier as a function of APC loss, f(1 − θP). Train on four regimes with θP at 1, 0.75, 0.5, 0.25, hold out three intermediate values, and the learned function reproduces the held-out truth to 0.12% with f(0) = 0 exactly and monotone throughout. Then freeze it and recover θP for a new regime. Two stages, calibrate then infer, and the function is learned across regimes instead of within one. The frozen inverse then gets 17, 15, 12, 9 of 35, about the same as the control, and θP in severe APC loss is still 127% off. A perfect learned function doesn't make its own input identifiable.
I'd summarize the hybrid results as: a network can learn a mechanism when you can vary the input to that mechanism across experiments, and when you can't, it learns whatever is cheapest, which is usually a constant that belongs to some other term.
depletion experiments
That summary suggested an experiment, so after the paper draft I ran it. If the network steals the basal constant because the data never gets near zero, add a condition that drives the regulator to zero and see if it gives the constant back.
Cheap version first: an equation-local screen that takes the exact reference states, fits only the equation with the learned term in it, holds everything outside that equation at truth, and starts the parameters inside it at 1.5× nominal. It strips out state estimation error, which makes it optimistic, but it runs in minutes instead of hours, so you can screen all twelve state-to-state edges across all four regimes. That atlas is where the "multiplicative modulators are easy" idea came from and then died: fine in normal, 45 to 65% functional error in the advanced and severe regimes for two of the four.
Then two depletion conditions. A retinoid-free condition takes RA to actual zero. A WNT-depletion condition lowers β-catenin but can't zero it, because HOXA13 feeds back into β-catenin production and holds a floor. RA → HOXA5 goes from 3 of 8 host-equation parameters recovered to 8 of 8 in the screen. RA → CYP26A1 goes 6 to 8. β-catenin → CYP26A1 gets worse, 30% to 65% basal error in severe, because β-catenin never reached zero and a_C is tiny so a small absolute offset is a big relative one. A twelve-condition control arm with two more ordinary conditions, same count, no depletion, moved nothing, so it's the direction of the intervention and not just having more data.
The sharper version is a dose sweep, condition count fixed at eleven, six doses, two regimes, two restarts. Three arms: RA restriction, WNT-only depletion, and WNT depletion plus a HOXA13 feedback knockdown that actually gets β-catenin to zero. At full dose the WNT-only arm leaves the β-catenin floor at about 0.01 of its max and a_M at 1.5 to 2% error. The combined arm reaches zero and a_M rounds to 0.0%. Same learned term, same seeds, the only difference is whether the anchor was actually reached.
I wrote down six predictions before running it and two passed: combined beats WNT-only at full dose, and the arms roughly agree at matched anchor ratios. Monotone in dose, an absolute error threshold for RA, and everything collapsing onto one curve did not hold. A five-start rerun of the severe half fixed the ugliest RA cell, 23% to 3.4%, and passed four of six, but RA at full depletion still lands at 3.3% basal error, which is exactly what the mechanistic control gets, so on that edge depletion isn't buying anything the control didn't already have.
Then the prospective test, which is where it got interesting. Take two hard edges, use the reference trajectories alone to design a protocol that reaches the anchor, run it. For MYC → HOXA13 the prescribed protocol drove MYC to zero and a_13 went from 12.8% error to 0.1% in normal and 95.6% to 0.0% in severe. The learned function got worse, 57% to 68% and 12% to 81%. The basal constant came back and the curve fell apart, partly because the input support changed under it. For HOXA13 → β-catenin, reaching the anchor required zeroing the WNT drive multiplier, which deletes W from the very condition that was supposed to identify W. W didn't move. So the rule is narrower than "observe near zero": you have to bring the learned function's input to zero while leaving the parameter you're trying to protect switched on, and for some edges in this model no condition does both. The reachability checker now rejects protocols that switch off the basal term. That's the actual finding from this section, and it's a negative one.
The full inverse PINN versions, twelve conditions with depletion, three starts, learned states and all: RA → HOXA5 gains eleven recovered parameters over its same-condition control across the four regimes, functional NRMSE 9.5 to 11.6%. β-catenin → MYC gets its curve down to 2.3 to 2.9% and its basal error stays at 10 to 17% because β-catenin still doesn't reach zero without the feedback knockdown, which that run didn't include. θP at high WNT is still wrong in every one of them.
All of this is one synthetic model with a known answer. The knockdowns are multipliers in a config file, not anything a lab has done.
compute
A 64-core box with one GPU, float64 throughout. Four regimes in parallel take 56 of the cores. There's no scheduler; the queue is a shell script in tmux that watches for a DONE line and launches the next variant. A hybrid variant, three restarts across four regimes, is about an hour and forty-five minutes. Inverse HMC is three and a half hours per regime. One edge atlas run died after 2.8 GPU-hours to a bug that would have been caught by a two-minute dry run, which I now do.
Things I got wrong along the way, for the record: the pseudoinverse Fisher matrix, the ESS-of-three posterior, the Sobol box that turned APC decay into growth, the stale forward plots, a cosine learning-rate schedule that froze the physical parameters before the networks had converged, and a best-of-n restart rule that picks by physics loss and in three regimes threw away the start with better parameter recovery. A third-stage refinement step with frozen networks is in every run and has never changed a result.
paper
The write-up is "A Reduced WNT–RA–HOX Model for Colorectal Cancer Stemness," with Nate and our advisor. It has the model, the nondimensionalization, the stemness score, and the three extensions above as sections. It was accepted to a SIAM conference in New York, and the project got a $6,000 internal research grant to keep going. Nate and I gave a 36-slide talk on it that took about 45 minutes, alternating. The depletion experiments happened after the draft and aren't in it yet.
The old ending of this post said PINNs are narrower than the enthusiasm suggests and good inside that range: known physics, not much data, a need to explain the answer. I'd keep that and add one thing. On an inverse problem, the method you use to fit the parameters matters much less than whether the parameters are identifiable at all, and PINNs make it very easy to not notice, because the trajectory always looks right. The most useful things I built this summer were the ones that tell you when not to trust the fit.
sources
Everything is in aidxhxr/PINN-Research: forward runs in PINN-smaller, the inverse lineage in PINN-inverse-*, the integral residual in PINN-inverse-pinn-boost, the Fisher and sensitivity scripts, inverse and forward HMC, and the hybrids and depletion experiments. Every number in this post has a run directory; this file maps them, including the ones I got wrong first. The figures are copied straight from those runs. The cover drawing and the slider above are illustrations.
The method is Raissi, Perdikaris, and Karniadakis. The Bayesian version follows Yang, Meng, and Karniadakis. The hybrids are universal differential equations by another name.