Can a Worm Play Chopin?

When you press a piano key:

How a neuron "decides":

Why C. elegans for music?

OpenWorm (2011–present) is an open-source project to recreate every cell of C. elegans in software:

# Step 2: K-means (Lab02 / L08 / L10 pattern)
from pyannow.step2_clustering.motor_primitives import (
    pca_reduce, find_motor_primitives, cluster_to_notes
)

scores, pca = pca_reduce(X_neural, n_components=4)
labels, km  = find_motor_primitives(scores, k=4)

# Each cluster transition → note onset
note_events = cluster_to_notes(labels, t_arr_ms)

# Step 3: Ridge (Lab07 / L07 / L11 pattern)
from pyannow.step3_regression.ridge_composer import RidgeComposer

ridge = RidgeComposer(alpha=None)  # CV chooses λ
ridge.fit(Z_worm, C_chopin)

C_pred   = ridge.predict(Z_worm)
eval_res = ridge.evaluate(Z_test, C_test)
print(f"Ridge R² = {eval_res['r2']:.3f}")

# Step 4-6: MLP + Adam + L-BFGS (Lab06/Lab10 pattern)
from pyannow.step4_ffnn.jax_composer import (
    create_model, init_params)
from pyannow.step5_training.adam_trainer import train_adam
from pyannow.step6_lbfgs.lbfgs_polish import polish_lbfgs

# Build MLP (same architecture as Lab10 MNIST)
model  = create_model(k_worm=4, k_chopin=8, hidden=32)
params = init_params(model, k_worm=4)

# Stage 1: Adam warmup (Lab10 / L20)
params, h_adam = train_adam(
    model, params, Z_worm, C_chopin,
    lr=1e-3, epochs=500, verbose=True
)

# Stage 2: L-BFGS polish (Lab08 / L22)
params, h_lbfgs = polish_lbfgs(
    model, params, Z_worm, C_chopin,
    max_steps=100
)
print(f"Final loss: {h_lbfgs.train_loss[-1]:.5f}")

Physics (ODE): each muscle group is a point oscillator

# Step 8a: ODE PINN (L14 + L27)
from pyannow.step8_pinn.locomotion_pinn import (
    PhysicsComposer, ode_physics_loss, train_pinn
)
import jax, jax.numpy as jnp

model_ode = PhysicsComposer(hidden=32, depth=2, out_dim=8)
params_ode = model_ode.init(
    jax.random.PRNGKey(0), jnp.zeros((1, 1+k_worm)))

# Collocation in TIME only (ODE)
t_colloc = jnp.linspace(0, 15, 200)
z_colloc = Z_worm[::T//200]

params_ode, hist_ode = train_pinn(
    model_ode, params_ode,
    inp_ode,       # (T, 1+k_worm)  ← time + neural code
    C_chopin,
    phys_loss_fn = ode_physics_loss,
    phys_args    = (t_colloc, z_colloc),
    lam_phys = 0.1,
    adam_steps=500, lbfgs_steps=50,
    label="ODE"
)

Physics (PDE): the worm body as a 1D elastic rod

# Step 8b: PDE PINN — 1D wave eq (L14 + L27)
from pyannow.step8_pinn.locomotion_pinn import (
    PhysicsComposer, pde_physics_loss,
    compare_ode_vs_pde
)

# Collocation in SPACE × TIME (PDE)
x_colloc = jnp.linspace(0, 1, 200)
t_colloc = jnp.linspace(0, 15, 200)

# One call runs BOTH and compares
results = compare_ode_vs_pde(
    Z_worm, C_chopin, t_arr,
    lam_phys=0.1,
    adam_steps=500, lbfgs_steps=50
)

print(f"ODE loss: {results['final_loss_ode']:.5f}")
print(f"PDE loss: {results['final_loss_pde']:.5f}")
# PDE is typically better — it knows the wave propagates!

The loss decomposition tells the story:

ODE physics loss: \(\|\ddot q + 2\gamma\dot q + \omega^2 q - F\|^2\) sampled at random times

PDE physics loss: \(\|\rho q_{tt} - \mu q_{xx} + \gamma q_t - F\|^2\) sampled at random (x, t) pairs

# Step 7: Random Forest (AppStat Lec07)
from sklearn.ensemble import RandomForestClassifier
from sklearn.preprocessing import StandardScaler

# Features: worm PCA (k) + beat position
X_rf = np.hstack([Z_worm_10s, beat_phase.reshape(-1,1)])
y_rf = onset_labels_10s  # binary onset labels

rf = RandomForestClassifier(
    n_estimators=100, max_depth=12,
    oob_score=True, random_state=0, n_jobs=-1)
rf.fit(X_rf_s, y_rf)

probs_rf  = rf.predict_proba(X_rf_s)[:,1]
onsets_rf = calibrated_onset_detect(probs_rf, ...)
F1_7      = musical_f1(onsets_rf, t_on_target)
print(f"Step 7 RF: F1={F1_7['f1']:.6f}")
print(f"  OOB score: {rf.oob_score_:.3f}")
# Feature importance: worm PC1 > beat_phase > PC2 …

# Step 9: Worm+Time MLP (NB04 recipe + AppStat MLP)
from sklearn.neural_network import MLPClassifier

def fourier_feats(t, T, n_harmonics=12):
    phase = t / T
    cols = [phase]
    for k in range(1, n_harmonics+1):
        cols += [np.cos(2*np.pi*k*phase),
                 np.sin(2*np.pi*k*phase)]
    return np.column_stack(cols)   # (T, 25-D)

T_feats = fourier_feats(t_arr, T=DURATION)   # 25-D
X9 = np.hstack([T_feats, Z_worm_scaled])      # 25+k

mlp9 = MLPClassifier(
    hidden_layer_sizes=(128,64), activation='relu',
    max_iter=300, early_stopping=True, random_state=0)
mlp9.fit(X9_sub, y9_sub, sample_weight=sw9)

F1_9 = musical_f1(onsets9, t_on_target)
# F1 = 0.933  (NB04 time-only ceiling: 0.858)
# Step 9b adds ODE residual → F1 ≈ 0.933 (tied)

Problem with Steps 9/9b: trained on 10s window, features extrapolate beyond training range for the full 229s piece → performance degrades.

# Step 10: NB04+Step9b full 229s piece — v1.0.0
t_on_all  = note_onsets(events_chopin, clip_s=10000)
T_full_s  = float(t_on_all[-1] + 2.0)   # ≈ 229s

# Run full Boyle model simulation
result_s10 = run_worm_simulation(T_full_s)
Z10 = result_s10['neural_pca']   # (T_pts, k_worm)

# Fourier features normalised over FULL piece
# → t/T_full ∈ [0,1] everywhere, never extrapolates
T_feats_full = fourier_feats(t_full, T=T_full_s)

# Stack: Fourier + worm PCA + ODE residual
X10 = np.hstack([T_feats_full, Z10, ode_res_full])

mlp10.fit(X10_sub, y10_sub, sample_weight=sw10)

# Synthesise and save full-length WAVs
audio_worm, fs = synthesise_melody(events10_full,
                                    duration_s=T_full_s)
wavfile.write('worm_mlp_full.wav', fs,
              (audio_worm * 28000).astype('int16'))

The gap today: PyANNOW uses synthetic neural activity (repeated muscle voltages). The real loop requires the actual Sibernetic SPH body simulator and C302 neural network running in a tight optimization cycle.

AMSC lectures directly applied: L03 (C++17), L05 (RAII subprocess), L09 (shared library), L11 (MPI reduce), L12 (OpenMP piano FEM), PC L5 (CUDA port of SPH).

The gap today: the PINN predicts (m_∞, τ_m, h_∞, τ_h) that minimise the PINN loss — but there is no statistical guarantee these predictions are biologically plausible relative to measured C. elegans channel values.

Full breakdown: AI_CONTRIBUTIONS.md