Quickstart Example - Membrane Dynamics with Nebuchadnezzar

Quickstart Example - Membrane Dynamics with Nebuchadnezzar

Overview

This quickstart demonstrates how to use the Membrane Dynamics package to create biophysically accurate circuit models for Nebuchadnezzar’s ATP-based differential equation system. We’ll model a simple membrane patch containing Na⁺/K⁺-ATPase pumps and voltage-gated ion channels, showing how membrane properties translate to circuit parameters and respond to ATP consumption.

Complete Working Example

use membrane_dynamics::prelude::*;
use nebuchadnezzar::prelude::*;

fn main() -> Result<(), Box<dyn std::error::Error>> {
    // 1. Create a membrane system with realistic composition
    let mut membrane = create_neuron_membrane_patch()?;
    
    // 2. Initialize ATP pool with physiological conditions
    let mut atp_pool = AtpPool::physiological();
    
    // 3. Set up membrane-circuit integration
    let mut circuit_interface = setup_membrane_circuit_interface()?;
    
    // 4. Run coupled membrane-circuit simulation
    let results = run_coupled_simulation(&mut membrane, &mut atp_pool, &mut circuit_interface)?;
    
    // 5. Analyze results
    analyze_and_display_results(&results);
    
    Ok(())
}

// Step 1: Create realistic neuron membrane patch
fn create_neuron_membrane_patch() -> Result<MembraneSystem, MembraneError> {
    // Define lipid composition (typical mammalian plasma membrane)
    let lipid_composition = LipidComposition::builder()
        .add_lipid(LipidType::PhosphatidylCholine, 0.45)    // 45% PC
        .add_lipid(LipidType::PhosphatidylSerine, 0.15)     // 15% PS  
        .add_lipid(LipidType::PhosphatidylEthanolamine, 0.20) // 20% PE
        .add_lipid(LipidType::Cholesterol, 0.20)            // 20% cholesterol
        .build();
    
    // Create membrane proteins with ATP dependencies
    let proteins = vec![
        // Na⁺/K⁺-ATPase (critical for membrane potential)
        MembraneProtein::builder()
            .protein_type(ProteinType::AtpPump {
                pump_type: PumpType::SodiumPotassium,
                stoichiometry: PumpStoichiometry {
                    atp_per_cycle: 1.0,
                    sodium_per_cycle: 3,   // 3 Na⁺ out
                    potassium_per_cycle: 2, // 2 K⁺ in
                },
                max_turnover_rate: 150.0, // s⁻¹
            })
            .density(100.0) // pumps/μm²
            .atp_km(0.5)    // mM - half-saturation for ATP
            .build(),
        
        // Voltage-gated sodium channels
        MembraneProtein::builder()
            .protein_type(ProteinType::IonChannel {
                ion_selectivity: IonSelectivity::Sodium,
                gating_mechanism: GatingMechanism::Voltage {
                    activation_voltage: -40.0, // mV
                    inactivation_voltage: -30.0, // mV
                    time_constant: 1.0,     // ms
                },
                conductance_range: (0.0, 20e-12), // 0-20 pS
            })
            .density(50.0)  // channels/μm²
            .build(),
        
        // Voltage-gated potassium channels
        MembraneProtein::builder()
            .protein_type(ProteinType::IonChannel {
                ion_selectivity: IonSelectivity::Potassium,
                gating_mechanism: GatingMechanism::Voltage {
                    activation_voltage: -50.0, // mV
                    inactivation_voltage: f64::INFINITY, // No inactivation
                    time_constant: 5.0,     // ms - slower than Na
                },
                conductance_range: (0.0, 15e-12), // 0-15 pS
            })
            .density(80.0)  // channels/μm²
            .build(),
    ];
    
    // Create membrane system
    let membrane = MembraneSystem::builder()
        .with_lipid_composition(lipid_composition)
        .with_proteins(proteins)
        .with_membrane_area(1e-12) // 1 μm² patch
        .with_temperature(310.0)   // 37°C
        .with_initial_voltage(-70.0) // -70 mV resting potential
        .build()?;
    
    Ok(membrane)
}

// Step 2: Set up membrane-circuit interface
fn setup_membrane_circuit_interface() -> Result<CircuitInterfaceLayer, InterfaceError> {
    let interface = CircuitInterfaceLayer::builder()
        // Enable all four membrane dynamics layers
        .with_molecular_layer(MolecularLayerConfig::full_physics())
        .with_mesoscale_layer(MesoscaleLayerConfig::with_domains())
        .with_cellular_layer(CellularLayerConfig::single_patch())
        .with_circuit_interface(CircuitInterfaceConfig::full_coupling())
        
        // Configure ATP coupling
        .with_atp_coupling_mode(AtpCouplingMode::Bidirectional)
        .with_atp_update_frequency(UpdateFrequency::EveryStep)
        
        // Set up hierarchical abstraction
        .with_hierarchical_mode(HierarchicalMode::Adaptive)
        .with_expansion_criteria(ExpansionCriteria {
            uncertainty_threshold: 0.1,    // Expand if uncertainty > 10%
            importance_threshold: 0.05,    // Expand if parameter importance > 5%
            computational_budget: 1000,    // Maximum circuit elements
        })
        
        .build()?;
    
    Ok(interface)
}

// Step 3: Run coupled simulation
fn run_coupled_simulation(
    membrane: &mut MembraneSystem,
    atp_pool: &mut AtpPool,
    circuit_interface: &mut CircuitInterfaceLayer,
) -> Result<SimulationResults, SimulationError> {
    // Initial membrane-to-circuit mapping
    let initial_circuit = circuit_interface.map_membrane_to_circuit(
        &membrane.state(),
        &atp_pool.state(),
    );
    
    // Initialize Nebuchadnezzar with membrane-derived circuit
    let mut nebuchadnezzar = NebuchadnezzarSystem::builder()
        .with_circuit_topology(initial_circuit)
        .with_atp_pool(atp_pool.clone())
        .with_solver(AtpSolver::AdaptiveRungeKutta {
            initial_dt: 0.001,  // 1 ms initial time step
            tolerance: 1e-6,    // Adaptive tolerance
        })
        .with_boundary_conditions(BoundaryConditions {
            external_voltage: 0.0,    // Ground reference
            ion_concentrations: IonConcentrations::physiological(),
        })
        .build()?;
    
    // Simulation parameters
    let total_time = 0.1;      // 100 ms simulation
    let recording_interval = 0.001; // Record every 1 ms
    let mut results = SimulationResults::new();
    
    // Main simulation loop
    let mut current_time = 0.0;
    while current_time < total_time {
        // Nebuchadnezzar ATP-based integration step
        let neb_result = nebuchadnezzar.integrate_step_adaptive()?;
        current_time += neb_result.actual_dt;
        
        // Extract ATP consumption information
        let atp_consumption = AtpConsumption {
            total_consumed: neb_result.atp_consumed,
            spatial_distribution: neb_result.spatial_atp_consumption,
            process_breakdown: neb_result.process_atp_breakdown,
        };
        
        // Update membrane based on ATP consumption
        let membrane_update = membrane.update_for_atp_consumption(
            &atp_consumption,
            neb_result.actual_dt,
        )?;
        
        // Update ATP pool
        atp_pool.consume_atp(atp_consumption.total_consumed);
        atp_pool.regenerate_atp(neb_result.actual_dt); // Background ATP synthesis
        
        // Update circuit parameters from membrane changes
        if membrane_update.significant_changes() {
            circuit_interface.update_circuit_parameters(
                &mut nebuchadnezzar.circuit,
                &membrane_update.parameter_changes,
            )?;
        }
        
        // Record results at specified intervals
        if (current_time / recording_interval).fract() < 0.01 { // Close to recording time
            results.record_time_point(TimePoint {
                time: current_time,
                membrane_state: membrane.state().clone(),
                atp_state: atp_pool.state().clone(),
                circuit_state: nebuchadnezzar.circuit.state().clone(),
                electrical_state: ElectricalState {
                    membrane_voltage: membrane.get_membrane_voltage(),
                    ion_currents: membrane.calculate_ion_currents(),
                    atp_consumption_rate: atp_consumption.total_consumed / neb_result.actual_dt,
                },
            });
        }
        
        // Adaptive hierarchy adjustment every 10 ms
        if (current_time * 100.0).round() % 10.0 == 0.0 {
            circuit_interface.adaptive_hierarchy_update(
                &nebuchadnezzar.get_sensitivity_analysis(),
            )?;
        }
    }
    
    Ok(results)
}

// Step 4: Analysis functions
fn analyze_and_display_results(results: &SimulationResults) {
    println!("=== Membrane Dynamics Simulation Results ===\n");
    
    // ATP efficiency analysis
    let atp_efficiency = calculate_atp_efficiency(results);
    println!("ATP Efficiency Metrics:");
    println!("  Average ATP consumption rate: {:.2} mM/s", atp_efficiency.avg_consumption_rate);
    println!("  Peak ATP consumption: {:.2} mM/s", atp_efficiency.peak_consumption);
    println!("  ATP utilization efficiency: {:.1}%", atp_efficiency.utilization_efficiency * 100.0);
    println!("  Energy cost per action potential: {:.3} fmol ATP\n", atp_efficiency.cost_per_action_potential);
    
    // Membrane parameter evolution
    let membrane_evolution = analyze_membrane_evolution(results);
    println!("Membrane Parameter Evolution:");
    println!("  Membrane capacitance range: {:.2} - {:.2} pF", 
             membrane_evolution.capacitance_range.0 * 1e12, 
             membrane_evolution.capacitance_range.1 * 1e12);
    println!("  Membrane resistance range: {:.1} - {:.1} GΩ", 
             membrane_evolution.resistance_range.0 / 1e9, 
             membrane_evolution.resistance_range.1 / 1e9);
    println!("  Voltage excursion: {:.1} to {:.1} mV", 
             membrane_evolution.voltage_range.0, 
             membrane_evolution.voltage_range.1);
    println!("  Na⁺/K⁺ pump activity range: {:.1} - {:.1} Hz\n", 
             membrane_evolution.pump_activity_range.0, 
             membrane_evolution.pump_activity_range.1);
    
    // Circuit parameter mapping validation
    let circuit_validation = validate_circuit_parameters(results);
    println!("Circuit Parameter Validation:");
    println!("  Membrane ↔ Circuit consistency: {:.1}%", circuit_validation.consistency_score * 100.0);
    println!("  ATP coupling accuracy: {:.1}%", circuit_validation.atp_coupling_accuracy * 100.0);
    println!("  Hierarchical abstraction error: {:.2}%", circuit_validation.abstraction_error * 100.0);
    println!("  Computational efficiency gain: {:.1}x\n", circuit_validation.efficiency_gain);
    
    // Biological realism assessment
    let realism_assessment = assess_biological_realism(results);
    println!("Biological Realism Assessment:");
    println!("  Resting potential accuracy: {:.1} mV (target: -70 mV)", realism_assessment.resting_potential);
    println!("  Action potential amplitude: {:.1} mV (target: ~110 mV)", realism_assessment.action_potential_amplitude);
    println!("  ATP consumption realism: {:.1}% of measured values", realism_assessment.atp_consumption_realism * 100.0);
    println!("  Ion pump stoichiometry accuracy: {:.1}%\n", realism_assessment.stoichiometry_accuracy * 100.0);
    
    // Generate plots if visualization is available
    #[cfg(feature = "plotting")]
    generate_analysis_plots(results);
}

// Analysis helper functions
fn calculate_atp_efficiency(results: &SimulationResults) -> AtpEfficiencyMetrics {
    let atp_consumption_rates: Vec<f64> = results.time_points
        .iter()
        .map(|tp| tp.electrical_state.atp_consumption_rate)
        .collect();
    
    AtpEfficiencyMetrics {
        avg_consumption_rate: atp_consumption_rates.iter().sum::<f64>() / atp_consumption_rates.len() as f64,
        peak_consumption: atp_consumption_rates.iter().fold(0.0, |a, &b| a.max(b)),
        utilization_efficiency: calculate_utilization_efficiency(results),
        cost_per_action_potential: estimate_action_potential_cost(results),
    }
}

fn analyze_membrane_evolution(results: &SimulationResults) -> MembraneEvolutionMetrics {
    let capacitances: Vec<f64> = results.time_points
        .iter()
        .map(|tp| tp.membrane_state.total_capacitance)
        .collect();
    
    let resistances: Vec<f64> = results.time_points
        .iter()
        .map(|tp| tp.membrane_state.total_resistance)
        .collect();
    
    let voltages: Vec<f64> = results.time_points
        .iter()
        .map(|tp| tp.electrical_state.membrane_voltage)
        .collect();
    
    MembraneEvolutionMetrics {
        capacitance_range: (
            capacitances.iter().fold(f64::INFINITY, |a, &b| a.min(b)),
            capacitances.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b)),
        ),
        resistance_range: (
            resistances.iter().fold(f64::INFINITY, |a, &b| a.min(b)),
            resistances.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b)),
        ),
        voltage_range: (
            voltages.iter().fold(f64::INFINITY, |a, &b| a.min(b)),
            voltages.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b)),
        ),
        pump_activity_range: calculate_pump_activity_range(results),
    }
}

// Optional plotting with visualization features
#[cfg(feature = "plotting")]
fn generate_analysis_plots(results: &SimulationResults) {
    use plotters::prelude::*;
    
    // Create output directory
    std::fs::create_dir_all("output/plots").unwrap();
    
    // Plot 1: Membrane voltage over time
    plot_membrane_voltage(results, "output/plots/membrane_voltage.png").unwrap();
    
    // Plot 2: ATP consumption rate over time
    plot_atp_consumption(results, "output/plots/atp_consumption.png").unwrap();
    
    // Plot 3: Circuit parameter evolution
    plot_circuit_parameters(results, "output/plots/circuit_evolution.png").unwrap();
    
    // Plot 4: Membrane-circuit correlation
    plot_membrane_circuit_correlation(results, "output/plots/correlation.png").unwrap();
    
    println!("Analysis plots saved to output/plots/");
}

// Data structures for results
#[derive(Debug, Clone)]
struct AtpEfficiencyMetrics {
    avg_consumption_rate: f64,
    peak_consumption: f64,
    utilization_efficiency: f64,
    cost_per_action_potential: f64,
}

#[derive(Debug, Clone)]
struct MembraneEvolutionMetrics {
    capacitance_range: (f64, f64),
    resistance_range: (f64, f64),
    voltage_range: (f64, f64),
    pump_activity_range: (f64, f64),
}

// Example command-line interface
#[cfg(feature = "cli")]
mod cli {
    use clap::{App, Arg};
    
    pub fn create_cli_app() -> App<'static, 'static> {
        App::new("Membrane Dynamics Example")
            .version("1.0")
            .author("Membrane Dynamics Team")
            .about("Demonstrates membrane-circuit coupling with Nebuchadnezzar")
            .arg(Arg::with_name("duration")
                .short("t")
                .long("time")
                .value_name("SECONDS")
                .help("Simulation duration in seconds")
                .takes_value(true)
                .default_value("0.1"))
            .arg(Arg::with_name("atp")
                .short("a")
                .long("atp-concentration")
                .value_name("MILLIMOLAR")
                .help("Initial ATP concentration in mM")
                .takes_value(true)
                .default_value("5.0"))
            .arg(Arg::with_name("output")
                .short("o")
                .long("output")
                .value_name("DIRECTORY")
                .help("Output directory for results")
                .takes_value(true)
                .default_value("output"))
            .arg(Arg::with_name("verbose")
                .short("v")
                .long("verbose")
                .help("Enable verbose output")
                .takes_value(false))
    }
}

Expected Output

When you run this example, you should see output similar to:

=== Membrane Dynamics Simulation Results ===

ATP Efficiency Metrics:
  Average ATP consumption rate: 0.42 mM/s
  Peak ATP consumption: 1.23 mM/s
  ATP utilization efficiency: 87.3%
  Energy cost per action potential: 2.15 fmol ATP

Membrane Parameter Evolution:
  Membrane capacitance range: 0.98 - 1.02 pF
  Membrane resistance range: 1.2 - 8.7 GΩ
  Voltage excursion: -75.3 to +45.2 mV
  Na⁺/K⁺ pump activity range: 45.2 - 189.7 Hz

Circuit Parameter Validation:
  Membrane ↔ Circuit consistency: 94.3%
  ATP coupling accuracy: 91.7%
  Hierarchical abstraction error: 2.1%
  Computational efficiency gain: 12.3x

Biological Realism Assessment:
  Resting potential accuracy: -69.8 mV (target: -70 mV)
  Action potential amplitude: 108.4 mV (target: ~110 mV)
  ATP consumption realism: 89.2% of measured values
  Ion pump stoichiometry accuracy: 97.8%

Analysis plots saved to output/plots/

Key Features Demonstrated

1. Realistic Membrane Composition: Physiologically accurate lipid and protein compositions
2. ATP-Dependent Dynamics: Membrane properties that change based on ATP consumption
3. Circuit Parameter Mapping: Direct translation from membrane biophysics to circuit elements
4. Bidirectional Coupling: ATP consumption affects membrane, membrane changes affect circuit
5. Hierarchical Abstraction: Adaptive computational optimization
6. Biological Validation: Results compared against experimental measurements

Next Steps

1. Experiment with different membrane compositions by modifying the create_neuron_membrane_patch() function
2. Add more complex proteins like ATP synthase or calcium pumps
3. Implement different cell types (cardiac, muscle, plant) with appropriate membrane properties
4. Explore optimization objectives beyond ATP efficiency
5. Integrate with experimental data for parameter validation and refinement

This example demonstrates the complete workflow for using Membrane Dynamics with Nebuchadnezzar to create biophysically accurate, ATP-driven circuit models of cellular membranes.

Membrane Dynamics Quickstart Example

This example demonstrates the complete membrane dynamics system integration with both the Nebuchadnezzar circuit system and the metacognitive orchestrator.

Complete Integration Example

1. Initialize Orchestrator-Managed Membrane System

import numpy as np
from membrane_dynamics import (
    MembranePatch, 
    CircuitInterface,
    OrchestratorInterface
)

# Initialize with orchestrator connection
orchestrator_interface = OrchestratorInterface(
    orchestrator_endpoint="ws://localhost:8888/orchestrator",
    module_id="membrane_dynamics_001"
)

# Create membrane patch under orchestrator supervision
membrane_patch = MembranePatch(
    area=1e-9,  # 1 μm² patch
    orchestrator=orchestrator_interface
)

# Circuit interface receives orchestrator directives
circuit_interface = CircuitInterface(
    nebuchadnezzar_endpoint="http://localhost:9999/circuits",
    orchestrator=orchestrator_interface
)

2. Orchestrator-Coordinated Membrane Setup

# Orchestrator provides system context
context_state = orchestrator_interface.get_current_context()
print(f"System ATP availability: {context_state['global_atp_pool']}")
print(f"Cognitive load priority: {context_state['processing_priority']}")

# Set up membrane components based on orchestrator guidance
atp_budget = orchestrator_interface.get_atp_allocation()

# Na⁺/K⁺-ATPase pumps (orchestrator manages ATP allocation)
na_k_atpase = ATPPump(
    density=1000,  # pumps/μm²
    max_atp_rate=atp_budget['na_k_pump_max'],  # orchestrator-limited
    orchestrator_managed=True
)

# Voltage-gated channels (orchestrator coordinates opening/closing)
vgsc = VoltageGatedChannel(
    type='Na', 
    density=500,
    orchestrator_prediction_enabled=True  # use intuition layer predictions
)

3. Real-time Orchestrator Coordination

def orchestrated_simulation_step(dt=0.01):
    """Single simulation step with orchestrator coordination"""
    
    # Receive orchestrator updates
    orchestrator_commands = orchestrator_interface.poll_commands()
    
    if 'context_update' in orchestrator_commands:
        membrane_patch.update_global_context(
            orchestrator_commands['context_update']
        )
    
    if 'reasoning_directive' in orchestrator_commands:
        circuit_interface.adjust_mapping_strategy(
            orchestrator_commands['reasoning_directive']
        )
    
    if 'intuition_prediction' in orchestrator_commands:
        membrane_patch.preload_predicted_changes(
            orchestrator_commands['intuition_prediction']
        )
    
    # Run membrane dynamics with orchestrator supervision
    membrane_state = membrane_patch.step(dt)
    
    # Report back to orchestrator
    status_report = {
        'membrane_voltage': membrane_state.voltage,
        'atp_consumption': membrane_state.atp_used,
        'circuit_parameters': circuit_interface.get_current_parameters(),
        'biological_realism_score': membrane_patch.validate_biology()
    }
    orchestrator_interface.report_status(status_report)
    
    return membrane_state

# Run orchestrated simulation
for t in np.arange(0, 10, 0.01):  # 10ms simulation
    state = orchestrated_simulation_step()
    
    # Orchestrator may adjust ATP allocation based on system needs
    if t > 5.0:  # Mid-simulation orchestrator intervention
        new_priority = orchestrator_interface.get_priority_update()
        if new_priority == 'cognitive_focus':
            # Reduce membrane ATP allocation for cognitive processing
            membrane_patch.reduce_maintenance_atp(factor=0.8)

4. Orchestrator-Validated Results