Biomedical

Protein folding under the Generalized Master Equation

Simulate non-Markovian folding dynamics with memory kernels. The Markov operator handles instantaneous transitions; the kernels carry the part of the past that still matters.

Source: examples/medicine_examples/protein_folding/main.rs

The full crate lives at examples/medicine_examples/protein_folding/. It is one of seven biomedical examples in the medicine crate, and the one that pairs cleanly with the Causal Monad: the simulation is a series of states, and each state depends on its predecessors in a way that ordinary Markov dynamics cannot capture.

What the model does

Protein folding is not memoryless. A protein that just passed through Intermediate 1 is more likely to keep moving toward the native state than a protein that arrived at Intermediate 1 by some other path. Ordinary Markov dynamics throw that information away. The Generalized Master Equation (GME) keeps it by adding a memory kernel:

P(t+Δt) = T · P(t)  +  Σ K(τ) · P(t−τ)

T is the instantaneous Markov transition matrix. K(τ) is the memory kernel at lag τ. P(t) is the probability distribution over conformational states. The kernel modifies the next step based on the recent history of the distribution.

Four conformational states

The example tracks four states: Unfolded, Intermediate 1, Intermediate 2, and Native (Folded). The initial distribution puts 100% of the probability mass at Unfolded. The Markov operator carries some of that mass forward each step; the memory kernels favor forward progression, so the path through I1 and I2 to Native gradually accumulates probability.

use deep_causality_core::EffectValue;
use deep_causality_physics::{Probability, generalized_master_equation};
use deep_causality_tensor::CausalTensor;

let num_states = 4;

let mut state: Vec<Probability> = vec![
    Probability::new(1.0).unwrap(),
    Probability::new(0.0).unwrap(),
    Probability::new(0.0).unwrap(),
    Probability::new(0.0).unwrap(),
];

main.rs:36. The Probability type from deep_causality_physics is a thin wrapper that constrains the inner f64 to [0, 1] at construction; later steps lean on that invariant.

The Markov operator

A 4×4 transition matrix encoded as a row-major Vec<f64>. Row i carries the probabilities of transitioning to state i from each of the four source states.

let markov_data = vec![
    // Row 0 (TO Unfolded): from U=0.7, from I1=0.1, from I2=0, from N=0
    0.70, 0.10, 0.00, 0.00,
    // Row 1 (TO I1): from U=0.3, from I1=0.7, from I2=0.1, from N=0
    0.30, 0.70, 0.10, 0.00,
    // Row 2 (TO I2): from U=0, from I1=0.2, from I2=0.4, from N=0
    0.00, 0.20, 0.40, 0.00,
    // Row 3 (TO Native): from U=0, from I1=0, from I2=0.5, from N=1.0
    0.00, 0.00, 0.50, 1.00,
];
let markov_operator = CausalTensor::new(markov_data, vec![num_states, num_states])?;

main.rs:55. The Native state is absorbing: row 3 has a 1.00 self-loop and no other inputs.

The memory kernels

Three kernels, one per lag. Each kernel is a small correction that nudges the next step toward the next forward conformation. The decay factor shrinks with lag, so the most recent past matters most.

let memory_kernels: Vec<CausalTensor<f64>> = (0..3)
    .map(|lag| {
        let decay = (-0.5 * (lag as f64 + 1.0)).exp() * 0.02;
        let mut data = vec![0.0; num_states * num_states];
        data[1] = decay;                 // Unfolded -> I1
        data[num_states + 2] = decay;    // I1 -> I2
        data[2 * num_states + 3] = decay; // I2 -> Native
        CausalTensor::new(data, vec![num_states, num_states]).unwrap()
    })
    .collect();

main.rs:71.

The simulation loop

Each tick applies generalized_master_equation to the current state, the history window, the Markov operator, and the kernels. The result is unwrapped, renormalized to handle floating-point drift, and pushed onto the history window.

for t in 1..=15 {
    let effect = generalized_master_equation(
        &state, &history, Some(&markov_operator), &memory_kernels,
    );
    if let EffectValue::Value(new_state) = effect.value() {
        // normalize, then update history (sliding window)
        history.remove(0);
        history.push(state.clone());
    }
}

The function returns a PropagatingEffect. Pattern-matching on EffectValue::Value gives you the next state vector; an error variant short-circuits without the surrounding loop having to know.

Run it

git clone https://github.com/deepcausality-rs/deep_causality
cd deep_causality
cargo run --release -p medicine_examples --example protein_folding

Expected output: a per-tick table showing how the probability mass migrates from Unfolded toward Native over 15 steps.

Why this is a good fit

Memory-dependent dynamics are the structural feature the Causal Monad was built to handle: each step is the same shape, the history is a typed value, and the composition produces a PropagatingEffect with the same audit-trail and error-handling guarantees as any other chain. Swapping in a different physical model is a function-pointer change, not a refactor.