Simulations
Four Simulations You Can Run Here
Each of these reproduces a number from the lesson. They run in your browser — Python compiled to WebAssembly, nothing installed, nothing sent anywhere. Press Run, change the constants, run again.
The figure above each runner is the expected result, so you can see what should happen before waiting for anything to load. The numbers to check: λ = 602 nm, a Landau–Zener agreement of a few hundredths, and a first-spike latency near 9 ms at 5 mW/mm².
⬇ Download the full script (numpy + scipy + matplotlib)The download uses scipy's adaptive integrator; the browser versions use a hand-written RK4, because scipy is not available in Pyodide. Simulation 2 says what that costs. Code is MIT-licensed.
Simulation 1 — Particle in a box
Steps 3–5 of the Part 1 derivation turned into code. For retinal (N = 12, L = 15.4 Å) it gives ΔE = 2.06 eV and λ = 602 nm, matching the hand calculation.
Try it: add one bond length to L to mimic wavefunction leakage at the ends, and see how far λ moves.
Simulation 2 — Landau–Zener crossing
The time-dependent Schrödinger equation for the two-state Hamiltonian of Part 2, with the diabatic energies swept linearly through each other (ℏ = 1).
Why this one disagrees slightly. The browser version integrates to a fixed time with a fixed step and lands within a few hundredths of the formula; the downloadable scipy version, with adaptive steps, reaches 0.014. The residual is not step size but the Stückelberg oscillation that has not yet died away at finite time — integrate further with the same step and it gets worse, because the states oscillate faster the further you go.
Try it: halve α (a slower wavepacket) and check that P falls as exp(−2πV²/α) predicts.
Simulation 3 — A light-driven neuron
A two-state model of channelrhodopsin (closed ⇆ open) coupled to the membrane equation of Part 3, with a leaky integrate-and-fire rule. The opening rate comes straight from the photon physics, kopen = σΦ × quantum yield.
\[ \frac{dp}{dt} = k_{\text{open}}(t)\,(1-p) - k_{\text{close}}\,p, \qquad C\,\frac{dV}{dt} = -\frac{V - V_{\text{rest}}}{R} + g_{\max}\,p\,(E_{\text{rev}} - V) \]
Try it: shorten the pulses to 1 ms, or raise the rate to 50 Hz. The channel's 10 ms closing time is the bottleneck.
Simulation 4 — Spike latency versus light intensity
The same model under continuous light, recording the time to the first spike across irradiances from 0.5 to 30 mW/mm².
Try it: overlay the analytical tth from Part 3, using I = gmax p∞ × 70 mV with p∞ = kopen / (kopen + kclose). Where and why does it differ?
Limits of these models
Real channelrhodopsin has at least four states, including a desensitised one that lowers the current during long pulses, and a real spike is shaped by voltage-gated Na⁺ and K⁺ channels rather than a reset rule. For research-grade work, swap in a four-state ChR2 model and a Hodgkin–Huxley neuron; the structure of the code stays the same.