diff --git a/build.rs b/build.rs index ad2e597..5b6218f 100644 --- a/build.rs +++ b/build.rs @@ -10,7 +10,8 @@ fn main() { println!("cargo:rerun-if-changed=src/ffi.rs"); println!("cargo:rerun-if-changed=cbindgen.toml"); - let crate_dir = PathBuf::from(env::var("CARGO_MANIFEST_DIR").expect("CARGO_MANIFEST_DIR not set")); + let crate_dir = + PathBuf::from(env::var("CARGO_MANIFEST_DIR").expect("CARGO_MANIFEST_DIR not set")); let crate_dir_string = crate_dir .to_str() .expect("crate directory must be valid UTF-8") @@ -36,7 +37,8 @@ fn main() { let mut generated = Vec::new(); bindings.write(&mut generated); - let header_contents = String::from_utf8(generated).expect("generated header was not valid UTF-8"); + let header_contents = + String::from_utf8(generated).expect("generated header was not valid UTF-8"); let header_path = header_dir.join("medicallib.h"); let needs_write = fs::read_to_string(&header_path) diff --git a/src/organs/bladder.rs b/src/organs/bladder.rs index b62a65b..e774d5f 100644 --- a/src/organs/bladder.rs +++ b/src/organs/bladder.rs @@ -1,21 +1,183 @@ use super::{Organ, OrganInfo}; use crate::types::OrganType; +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum BladderPhase { + Filling, + Voiding, + PostVoidRefractory, +} + #[derive(Debug, Clone)] pub struct Bladder { info: OrganInfo, - /// Urine volume ml + /// Urine volume stored in the bladder (ml). pub volume_ml: f32, - /// Pressure proxy (cmH2O) + /// Intraluminal/detrusor pressure (cmH2O). pub pressure: f32, + /// Normalized afferent firing (0..=1) representing stretch receptor activity. + pub afferent_signal: f32, + /// Normalized urgency perception (0..=1). + pub urgency: f32, + /// Current phase of the bladder state machine. + pub phase: BladderPhase, + /// Functional capacity where continence is expected (ml). + pub capacity_ml: f32, + /// Residual volume expected after a complete void (ml). + pub residual_volume_ml: f32, + /// Compliance relating volume change to pressure (ml per cmH2O). + pub compliance_ml_per_cm_h2o: f32, + /// Baseline pressure generated by abdominal cavity (cmH2O). + pub baseline_pressure_cm_h2o: f32, + /// Volume threshold at which a full micturition reflex is triggered (ml). + pub micturition_threshold_ml: f32, + /// Volume threshold where urge perception begins (ml). + pub urge_threshold_ml: f32, + /// Average renal inflow into the bladder (ml/min). + pub filling_rate_ml_per_min: f32, + /// Peak voluntary/automatic outflow during voiding (ml/s). + pub voiding_flow_ml_per_s: f32, + /// Tone of the internal urethral sphincter (0..=1, higher means more closed). + pub internal_sphincter_tone: f32, + /// Tone of the external urethral sphincter/pelvic floor (0..=1). + pub external_sphincter_tone: f32, + /// Parasympathetic drive to the detrusor (0..=1). + pub parasympathetic_drive: f32, + /// Sympathetic drive maintaining storage (0..=1). + pub sympathetic_drive: f32, + /// Somatic drive through the pudendal nerve to the external sphincter (0..=1). + pub somatic_drive: f32, + time_since_last_void_s: f32, + refractory_seconds: f32, } impl Bladder { pub fn new(id: impl Into) -> Self { Self { info: OrganInfo::new(id, OrganType::Bladder), - volume_ml: 100.0, + volume_ml: 120.0, pressure: 5.0, + afferent_signal: 0.1, + urgency: 0.1, + phase: BladderPhase::Filling, + capacity_ml: 500.0, + residual_volume_ml: 30.0, + compliance_ml_per_cm_h2o: 18.0, + baseline_pressure_cm_h2o: 5.0, + micturition_threshold_ml: 350.0, + urge_threshold_ml: 200.0, + filling_rate_ml_per_min: 60.0, + voiding_flow_ml_per_s: 15.0, + internal_sphincter_tone: 0.85, + external_sphincter_tone: 0.9, + parasympathetic_drive: 0.05, + sympathetic_drive: 0.8, + somatic_drive: 0.8, + time_since_last_void_s: 0.0, + refractory_seconds: 15.0, + } + } + + fn approach(current: f32, target: f32, rate_per_second: f32, dt_seconds: f32) -> f32 { + let rate = rate_per_second.max(0.0); + if rate == 0.0 || dt_seconds <= 0.0 { + return current; + } + let delta = target - current; + let max_step = rate * dt_seconds; + if delta > max_step { + current + max_step + } else if delta < -max_step { + current - max_step + } else { + target + } + } + + fn update_drives(&mut self, dt_seconds: f32) { + let urgency = self.urgency; + let (parasym_target, sym_target, somatic_target) = match self.phase { + BladderPhase::Filling => { + let parasym = urgency.powf(1.2).clamp(0.0, 0.85); + let sym = (1.0 - 0.6 * urgency).clamp(0.2, 1.0); + let somatic = (0.95 - 0.5 * urgency).clamp(0.3, 0.95); + (parasym, sym, somatic) + } + BladderPhase::Voiding => (1.0, 0.1, 0.2), + BladderPhase::PostVoidRefractory => (0.2, 0.6, 0.7), + }; + + self.parasympathetic_drive = + Self::approach(self.parasympathetic_drive, parasym_target, 0.8, dt_seconds); + self.sympathetic_drive = + Self::approach(self.sympathetic_drive, sym_target, 0.6, dt_seconds); + self.somatic_drive = Self::approach(self.somatic_drive, somatic_target, 0.9, dt_seconds); + + self.parasympathetic_drive = self.parasympathetic_drive.clamp(0.0, 1.0); + self.sympathetic_drive = self.sympathetic_drive.clamp(0.0, 1.0); + self.somatic_drive = self.somatic_drive.clamp(0.0, 1.0); + } + + fn update_sphincters(&mut self) { + let internal = 0.3 + 0.7 * self.sympathetic_drive; + let external = 0.2 + 0.8 * self.somatic_drive; + self.internal_sphincter_tone = internal.clamp(0.0, 1.0); + self.external_sphincter_tone = external.clamp(0.0, 1.0); + } + + fn update_afferents(&mut self) { + let low_volume_component = (self.volume_ml / 50.0).clamp(0.0, 1.0) * 0.15; + let fullness_component = if self.capacity_ml > self.urge_threshold_ml { + let denom = (self.capacity_ml - self.urge_threshold_ml).max(1.0); + ((self.volume_ml - self.urge_threshold_ml) / denom).clamp(0.0, 1.0) + } else { + (self.volume_ml / self.capacity_ml.max(1.0)).clamp(0.0, 1.0) + }; + self.afferent_signal = (low_volume_component + fullness_component).clamp(0.0, 1.0); + self.urgency = self.afferent_signal.powf(1.35).clamp(0.0, 1.0); + } + + fn update_pressure(&mut self) { + let passive_volume_ml = (self.volume_ml - 30.0).max(0.0); + let normalized_volume = (self.volume_ml / self.capacity_ml).clamp(0.0, 1.5); + let compliance_factor = 1.0 + 4.0 * normalized_volume.powf(4.0); + let passive_pressure = if self.compliance_ml_per_cm_h2o > 0.0 { + passive_volume_ml / self.compliance_ml_per_cm_h2o * compliance_factor + } else { + 0.0 + }; + let active_pressure = 40.0 * self.parasympathetic_drive; + let abdominal = self.baseline_pressure_cm_h2o; + self.pressure = (abdominal + passive_pressure + active_pressure).clamp(0.0, 90.0); + } + + fn handle_filling_phase(&mut self, _dt_seconds: f32) { + if self.volume_ml >= self.micturition_threshold_ml || self.pressure > 45.0 { + if self.external_sphincter_tone < 0.4 || self.urgency > 0.95 { + self.phase = BladderPhase::Voiding; + } + } + let overdistention_limit = self.capacity_ml * 1.4; + if self.volume_ml > overdistention_limit { + self.volume_ml = overdistention_limit; + } + } + + fn handle_voiding_phase(&mut self, dt_seconds: f32) { + let relaxation_factor = 1.0 - 0.5 * self.external_sphincter_tone; + let drive = (self.parasympathetic_drive * relaxation_factor).clamp(0.0, 1.0); + let pressure_factor = (self.pressure / 40.0).clamp(0.0, 2.5); + let outflow = self.voiding_flow_ml_per_s.max(0.0) * pressure_factor * drive * dt_seconds; + self.volume_ml = (self.volume_ml - outflow).max(self.residual_volume_ml); + if self.volume_ml <= self.residual_volume_ml + 1.0 { + self.phase = BladderPhase::PostVoidRefractory; + self.time_since_last_void_s = 0.0; + } + } + + fn handle_post_void_phase(&mut self) { + if self.time_since_last_void_s >= self.refractory_seconds { + self.phase = BladderPhase::Filling; } } } @@ -27,16 +189,42 @@ impl Organ for Bladder { fn organ_type(&self) -> OrganType { self.info.kind() } - fn update(&mut self, _dt_seconds: f32) { - // simplistic pressure-volume relation - self.pressure = (self.volume_ml / 50.0).clamp(0.0, 30.0); + fn update(&mut self, dt_seconds: f32) { + if dt_seconds <= 0.0 { + return; + } + + self.time_since_last_void_s += dt_seconds; + let inflow = (self.filling_rate_ml_per_min / 60.0).max(0.0) * dt_seconds; + self.volume_ml += inflow; + + self.update_afferents(); + self.update_drives(dt_seconds); + self.update_sphincters(); + self.update_pressure(); + + match self.phase { + BladderPhase::Filling => self.handle_filling_phase(dt_seconds), + BladderPhase::Voiding => self.handle_voiding_phase(dt_seconds), + BladderPhase::PostVoidRefractory => self.handle_post_void_phase(), + } + + if matches!(self.phase, BladderPhase::Voiding) + && self.volume_ml <= self.residual_volume_ml + 1.0 + { + self.phase = BladderPhase::PostVoidRefractory; + self.time_since_last_void_s = 0.0; + } } fn summary(&self) -> String { format!( - "Bladder[id={}, vol={:.0} ml, P={:.1}]", + "Bladder[id={}, phase={:?}, vol={:.0}/{:.0} ml, P={:.1} cmH2O, urge={:.0}%]", self.id(), + self.phase, self.volume_ml, - self.pressure + self.capacity_ml, + self.pressure, + self.urgency * 100.0 ) } fn as_any(&self) -> &dyn core::any::Any { diff --git a/src/organs/brain.rs b/src/organs/brain.rs index 78537d5..76ffc7e 100644 --- a/src/organs/brain.rs +++ b/src/organs/brain.rs @@ -1,21 +1,173 @@ use super::{Organ, OrganInfo}; use crate::types::OrganType; +use core::f32::consts::TAU; + +const CIRCADIAN_PERIOD_SECONDS: f32 = 24.0 * 3600.0; +const HOMEOSTATIC_WAKE_ACCUMULATION_S: f32 = 16.0 * 3600.0; +const HOMEOSTATIC_DISCHARGE_S: f32 = 6.0 * 3600.0; + +/// Sleep architecture stages used for brain state transitions. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum SleepStage { + Wake, + N1, + N2, + N3, + Rem, +} #[derive(Debug, Clone)] pub struct Brain { info: OrganInfo, - /// 0..=100 scale of consciousness. + /// 0..=100 index approximating global cortical consciousness/alertness. pub consciousness: u8, - /// Simplified neural activity index. + /// Composite neural activity index blending frequency and metabolic power. pub activity_index: f32, + /// Normalized cortical arousal (0..=1). + pub cortical_arousal: f32, + /// Brainstem autonomic output regulating MAP and respiratory drive (0..=1). + pub brainstem_autonomic_drive: f32, + /// Intracranial pressure (mmHg). + pub intracranial_pressure_mm_hg: f32, + /// Cerebral perfusion pressure (mmHg). + pub cerebral_perfusion_pressure_mm_hg: f32, + /// Cerebral blood flow (ml/100g/min). + pub cerebral_blood_flow_ml_per_100g_min: f32, + /// Fractional metabolic demand relative to resting wakefulness. + pub metabolic_demand_fraction: f32, + /// Cortical oxygen saturation (fraction 0..=1). + pub oxygenation_saturation: f32, + /// Excitatory glutamatergic tone (relative units 0.3..=1.4). + pub glutamate_level: f32, + /// Inhibitory GABAergic tone (relative units 0.4..=1.2). + pub gaba_level: f32, + /// Mesolimbic/striatal dopamine tone (0..=1). + pub dopamine_tone: f32, + /// Sleep homeostatic pressure (0..≈1.1). + pub sleep_pressure: f32, + /// Circadian phase in radians (0..TAU, midnight ≈ 0). + pub circadian_phase_radians: f32, + /// Current polysomnographic sleep stage. + pub sleep_stage: SleepStage, + /// Dominant EEG frequency (Hz). + pub eeg_dominant_frequency_hz: f32, + /// Instantaneous seizure risk (0..=1). + pub seizure_risk: f32, + /// Autonomic variability (0..=1, higher reflects sympathetic swings). + pub autonomic_variability: f32, + /// Cognitive/task load (0..=1) influencing metabolic demand. + pub cognitive_load: f32, + time_in_stage_s: f32, } impl Brain { pub fn new(id: impl Into) -> Self { Self { info: OrganInfo::new(id, OrganType::Brain), - consciousness: 100, + consciousness: 95, activity_index: 1.0, + cortical_arousal: 0.82, + brainstem_autonomic_drive: 0.9, + intracranial_pressure_mm_hg: 10.0, + cerebral_perfusion_pressure_mm_hg: 75.0, + cerebral_blood_flow_ml_per_100g_min: 52.0, + metabolic_demand_fraction: 1.05, + oxygenation_saturation: 0.98, + glutamate_level: 0.65, + gaba_level: 0.55, + dopamine_tone: 0.6, + sleep_pressure: 0.4, + circadian_phase_radians: TAU * 0.25, + sleep_stage: SleepStage::Wake, + eeg_dominant_frequency_hz: 18.0, + seizure_risk: 0.05, + autonomic_variability: 0.35, + cognitive_load: 0.35, + time_in_stage_s: 0.0, + } + } + + fn approach(current: f32, target: f32, rate_per_second: f32, dt_seconds: f32) -> f32 { + let rate = rate_per_second.max(0.0); + if rate == 0.0 || dt_seconds <= 0.0 { + return current; + } + let delta = target - current; + let max_step = rate * dt_seconds; + if delta > max_step { + current + max_step + } else if delta < -max_step { + current - max_step + } else { + target + } + } + + fn wrap_phase(phase: f32) -> f32 { + if phase >= 0.0 && phase < TAU { + return phase; + } + let mut wrapped = phase % TAU; + if wrapped < 0.0 { + wrapped += TAU; + } + wrapped + } + + fn transition_stage(&mut self, stage: SleepStage) { + if self.sleep_stage != stage { + self.sleep_stage = stage; + self.time_in_stage_s = 0.0; + } + } + + fn evaluate_sleep_stage(&mut self, base_arousal: f32) { + match self.sleep_stage { + SleepStage::Wake => { + if base_arousal < 0.35 && self.sleep_pressure > 0.55 { + self.transition_stage(SleepStage::N1); + } + } + SleepStage::N1 => { + if base_arousal > 0.5 { + self.transition_stage(SleepStage::Wake); + } else if self.time_in_stage_s > 180.0 && self.sleep_pressure > 0.6 { + self.transition_stage(SleepStage::N2); + } + } + SleepStage::N2 => { + if base_arousal > 0.52 { + self.transition_stage(SleepStage::Wake); + } else if self.time_in_stage_s > 900.0 && self.sleep_pressure > 0.65 { + self.transition_stage(SleepStage::N3); + } else if self.time_in_stage_s > 2400.0 { + self.transition_stage(SleepStage::Rem); + } + } + SleepStage::N3 => { + if self.time_in_stage_s > 1800.0 { + self.transition_stage(SleepStage::Rem); + } else if base_arousal > 0.45 { + self.transition_stage(SleepStage::N2); + } + } + SleepStage::Rem => { + if self.time_in_stage_s > 1500.0 { + self.transition_stage(SleepStage::N2); + } else if base_arousal > 0.65 { + self.transition_stage(SleepStage::Wake); + } + } + } + } + + fn stage_arousal_target(&self, base_arousal: f32) -> f32 { + match self.sleep_stage { + SleepStage::Wake => base_arousal.max(0.6), + SleepStage::N1 => base_arousal.clamp(0.3, 0.55), + SleepStage::N2 => base_arousal.clamp(0.2, 0.45), + SleepStage::N3 => base_arousal.min(0.25), + SleepStage::Rem => base_arousal.clamp(0.45, 0.7), } } } @@ -27,15 +179,195 @@ impl Organ for Brain { fn organ_type(&self) -> OrganType { self.info.kind() } - fn update(&mut self, _dt_seconds: f32) { - self.activity_index = self.activity_index.clamp(0.0, 2.0); + fn update(&mut self, dt_seconds: f32) { + if dt_seconds <= 0.0 { + return; + } + + self.time_in_stage_s += dt_seconds; + let circadian_increment = TAU * (dt_seconds / CIRCADIAN_PERIOD_SECONDS); + self.circadian_phase_radians = + Self::wrap_phase(self.circadian_phase_radians + circadian_increment); + + let circadian_arousal = 0.55 + 0.45 * (self.circadian_phase_radians - TAU * 0.25).sin(); + + let asleep = !matches!(self.sleep_stage, SleepStage::Wake); + let sleep_pressure_delta = if asleep { + -dt_seconds / HOMEOSTATIC_DISCHARGE_S + } else { + dt_seconds / HOMEOSTATIC_WAKE_ACCUMULATION_S + }; + self.sleep_pressure = (self.sleep_pressure + sleep_pressure_delta).clamp(0.0, 1.1); + + let base_arousal = (circadian_arousal - 0.55 * self.sleep_pressure).clamp(0.05, 1.0); + self.evaluate_sleep_stage(base_arousal); + + let arousal_target = self.stage_arousal_target(base_arousal); + self.cortical_arousal = + Self::approach(self.cortical_arousal, arousal_target, 0.8, dt_seconds).clamp(0.05, 1.0); + + let (autonomic_target, variability_target) = match self.sleep_stage { + SleepStage::Wake => (0.9, 0.35), + SleepStage::N1 => (0.82, 0.4), + SleepStage::N2 => (0.78, 0.45), + SleepStage::N3 => (0.72, 0.3), + SleepStage::Rem => (0.86, 0.6), + }; + self.brainstem_autonomic_drive = Self::approach( + self.brainstem_autonomic_drive, + autonomic_target, + 0.6, + dt_seconds, + ) + .clamp(0.3, 1.1); + self.autonomic_variability = Self::approach( + self.autonomic_variability, + variability_target, + 0.4, + dt_seconds, + ) + .clamp(0.05, 0.95); + + let cognitive_target = if matches!(self.sleep_stage, SleepStage::Wake) { + (0.3 + 0.45 * self.cortical_arousal).clamp(0.15, 1.0) + } else { + 0.12 + }; + self.cognitive_load = + Self::approach(self.cognitive_load, cognitive_target, 0.25, dt_seconds) + .clamp(0.05, 1.0); + + let base_metabolic = match self.sleep_stage { + SleepStage::Wake => 1.05, + SleepStage::N1 => 0.95, + SleepStage::N2 => 0.85, + SleepStage::N3 => 0.7, + SleepStage::Rem => 1.0, + }; + let metabolic_target = (base_metabolic + + 0.25 * (self.cognitive_load - 0.2) + + 0.1 * (self.glutamate_level - self.gaba_level)) + .clamp(0.6, 1.4); + self.metabolic_demand_fraction = Self::approach( + self.metabolic_demand_fraction, + metabolic_target, + 0.35, + dt_seconds, + ) + .clamp(0.6, 1.5); + + let (glu_base, gaba_base, dopamine_base, eeg_target_base) = match self.sleep_stage { + SleepStage::Wake => (0.65, 0.55, 0.6, 18.0), + SleepStage::N1 => (0.55, 0.62, 0.5, 9.0), + SleepStage::N2 => (0.5, 0.68, 0.45, 6.0), + SleepStage::N3 => (0.45, 0.75, 0.4, 2.0), + SleepStage::Rem => (0.68, 0.6, 0.65, 10.5), + }; + let arousal_offset = self.cortical_arousal - 0.5; + let glutamate_target = (glu_base + 0.15 * arousal_offset).clamp(0.3, 1.3); + let gaba_target = (gaba_base - 0.1 * arousal_offset).clamp(0.4, 1.2); + let dopamine_target = (dopamine_base + 0.05 * circadian_arousal + - 0.03 * self.sleep_pressure) + .clamp(0.35, 0.85); + let eeg_target = match self.sleep_stage { + SleepStage::Wake => eeg_target_base + 4.0 * (self.cortical_arousal - 0.7).max(0.0), + SleepStage::Rem => eeg_target_base + 2.0 * (self.cortical_arousal - 0.6).max(0.0), + _ => eeg_target_base, + }; + + self.glutamate_level = + Self::approach(self.glutamate_level, glutamate_target, 0.5, dt_seconds).clamp(0.3, 1.4); + self.gaba_level = + Self::approach(self.gaba_level, gaba_target, 0.45, dt_seconds).clamp(0.4, 1.2); + self.dopamine_tone = + Self::approach(self.dopamine_tone, dopamine_target, 0.3, dt_seconds).clamp(0.3, 0.9); + self.eeg_dominant_frequency_hz = + Self::approach(self.eeg_dominant_frequency_hz, eeg_target, 0.8, dt_seconds) + .clamp(0.5, 35.0); + + let oxygen_target = (0.95 + 0.03 * self.brainstem_autonomic_drive + - 0.04 * (self.metabolic_demand_fraction - 1.0)) + .clamp(0.88, 0.99); + self.oxygenation_saturation = + Self::approach(self.oxygenation_saturation, oxygen_target, 0.4, dt_seconds) + .clamp(0.8, 1.0); + + let cbf_target = (50.0 * self.metabolic_demand_fraction + + 0.25 * (self.cerebral_perfusion_pressure_mm_hg - 70.0)) + .clamp(30.0, 90.0); + self.cerebral_blood_flow_ml_per_100g_min = Self::approach( + self.cerebral_blood_flow_ml_per_100g_min, + cbf_target, + 1.6, + dt_seconds, + ) + .clamp(25.0, 95.0); + + let stage_icp_term = match self.sleep_stage { + SleepStage::Wake => 0.0, + SleepStage::N1 => 0.5, + SleepStage::N2 => 1.0, + SleepStage::N3 => 1.8, + SleepStage::Rem => 1.2, + }; + let icp_target = (10.0 + + stage_icp_term + + 0.12 * (self.cerebral_blood_flow_ml_per_100g_min - 50.0) + + 4.0 * (self.metabolic_demand_fraction - 1.0)) + .clamp(5.0, 30.0); + self.intracranial_pressure_mm_hg = Self::approach( + self.intracranial_pressure_mm_hg, + icp_target, + 0.4, + dt_seconds, + ) + .clamp(4.0, 35.0); + + let map_proxy = 90.0 + 15.0 * (self.brainstem_autonomic_drive - 0.8) + - 8.0 * (self.autonomic_variability - 0.4); + let cpp_target = (map_proxy - self.intracranial_pressure_mm_hg).clamp(40.0, 110.0); + self.cerebral_perfusion_pressure_mm_hg = Self::approach( + self.cerebral_perfusion_pressure_mm_hg, + cpp_target, + 1.2, + dt_seconds, + ); + + let activity = (self.cortical_arousal * 0.6 + + self.metabolic_demand_fraction * 0.25 + + (self.glutamate_level - self.gaba_level + 0.7) * 0.1 + + self.dopamine_tone * 0.05) + .clamp(0.1, 2.3); + self.activity_index = activity; + + let perfusion_factor = (self.cerebral_perfusion_pressure_mm_hg / 70.0).clamp(0.4, 1.3); + let oxygen_factor = (self.oxygenation_saturation / 0.96).clamp(0.5, 1.1); + let stage_bonus = match self.sleep_stage { + SleepStage::Wake => 18.0, + SleepStage::Rem => 8.0, + _ => 3.0, + }; + let target_consciousness = + ((self.cortical_arousal * 70.0 * perfusion_factor * oxygen_factor) + stage_bonus) + .clamp(0.0, 100.0); + self.consciousness = target_consciousness.round() as u8; + + let excitability = + (self.glutamate_level - self.gaba_level + self.cortical_arousal - 0.5).max(0.0); + let hypoxia = (0.94 - self.oxygenation_saturation).max(0.0) * 4.0; + let perfusion_deficit = (55.0 - self.cerebral_perfusion_pressure_mm_hg).max(0.0) / 40.0; + self.seizure_risk = (0.25 * excitability + hypoxia + perfusion_deficit).clamp(0.0, 1.0); } fn summary(&self) -> String { format!( - "Brain[id={}, GCS~{}, activity={:.2}]", + "Brain[id={}, stage={:?}, arousal={:.2}, ICP={:.1} mmHg, CPP={:.0} mmHg, O2={:.0}%, szrisk={:.0}%]", self.id(), - self.consciousness, - self.activity_index + self.sleep_stage, + self.cortical_arousal, + self.intracranial_pressure_mm_hg, + self.cerebral_perfusion_pressure_mm_hg, + self.oxygenation_saturation * 100.0, + self.seizure_risk * 100.0 ) } fn as_any(&self) -> &dyn core::any::Any { diff --git a/src/organs/esophagus.rs b/src/organs/esophagus.rs index 9c557f6..020044a 100644 --- a/src/organs/esophagus.rs +++ b/src/organs/esophagus.rs @@ -1,20 +1,253 @@ use super::{Organ, OrganInfo}; use crate::types::OrganType; +const ESOPHAGEAL_LENGTH_CM: f32 = 25.0; +const PRIMARY_WAVE_SPEED_CM_S: f32 = 4.5; +const SECONDARY_WAVE_SPEED_CM_S: f32 = 3.2; +const BASE_HIATAL_PRESSURE_CM_H2O: f32 = 6.0; + +/// Functional sequence of the esophagus during swallowing and reflux handling. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum EsophagealStage { + Idle, + SwallowInitiation, + PrimaryPeristalsis, + SecondaryPeristalsis, + Clearing, + RefluxExposure, +} + #[derive(Debug, Clone)] pub struct Esophagus { info: OrganInfo, - /// Reflux severity 0..=100 - pub reflux: u8, + /// Luminal pH along the distal esophagus. + pub luminal_ph: f32, + /// Lower esophageal sphincter tone (0..=1). + pub lower_sphincter_tone: f32, + /// Upper esophageal sphincter tone (0..=1). + pub upper_sphincter_tone: f32, + /// Distance traversed by the current peristaltic wave (cm). + pub peristaltic_progress_cm: f32, + /// Bolus volume currently descending (ml). + pub bolus_volume_ml: f32, + /// Salivary buffer content available to neutralize acid (ml). + pub saliva_buffer_ml: f32, + /// Estimated reflux episodes per hour. + pub reflux_events_per_hour: f32, + /// Fractional acid exposure burden (0..≈1.2). + pub acid_exposure_fraction: f32, + /// Mucosal integrity (0..≈1; <0.7 suggests erosive disease). + pub mucosal_integrity: f32, + /// Current functional state. + pub stage: EsophagealStage, + /// Swallow drive (0..=1) influenced by salivary demand and mucosal irritation. + pub swallow_drive: f32, + /// Peristaltic contractile strength scaling factor. + pub peristaltic_strength: f32, + /// Vagal tone modulating motility (0..=1). + pub vagal_tone: f32, + /// Time since the last initiated swallow (s). + pub time_since_last_swallow_s: f32, + time_in_stage_s: f32, + /// Target interval between swallows based on drive (s). + pub swallow_interval_target_s: f32, + /// Pressure gradient promoting reflux (cmH2O). + pub hiatal_pressure_gradient_cm_h2o: f32, } impl Esophagus { pub fn new(id: impl Into) -> Self { Self { info: OrganInfo::new(id, OrganType::Esophagus), - reflux: 0, + luminal_ph: 6.4, + lower_sphincter_tone: 0.78, + upper_sphincter_tone: 0.9, + peristaltic_progress_cm: 0.0, + bolus_volume_ml: 0.0, + saliva_buffer_ml: 1.2, + reflux_events_per_hour: 1.5, + acid_exposure_fraction: 0.08, + mucosal_integrity: 0.95, + stage: EsophagealStage::Idle, + swallow_drive: 0.35, + peristaltic_strength: 0.9, + vagal_tone: 0.65, + time_since_last_swallow_s: 0.0, + time_in_stage_s: 0.0, + swallow_interval_target_s: 22.0, + hiatal_pressure_gradient_cm_h2o: BASE_HIATAL_PRESSURE_CM_H2O, } } + + fn approach(current: f32, target: f32, rate_per_second: f32, dt_seconds: f32) -> f32 { + let rate = rate_per_second.max(0.0); + if rate == 0.0 || dt_seconds <= 0.0 { + return current; + } + let delta = target - current; + let max_step = rate * dt_seconds; + if delta > max_step { + current + max_step + } else if delta < -max_step { + current - max_step + } else { + target + } + } + + fn transition_stage(&mut self, stage: EsophagealStage) { + if self.stage != stage { + self.stage = stage; + self.time_in_stage_s = 0.0; + } + } + + fn start_primary_swallow(&mut self) { + self.transition_stage(EsophagealStage::SwallowInitiation); + self.peristaltic_progress_cm = 0.0; + self.bolus_volume_ml = 4.0 + 3.0 * self.swallow_drive; + self.saliva_buffer_ml += 1.5 + 1.2 * self.swallow_drive; + self.time_since_last_swallow_s = 0.0; + } + + fn update_swallow_drive(&mut self, dt_seconds: f32) { + let dryness = (self.time_since_last_swallow_s / 45.0).clamp(0.0, 1.4); + let irritation = (self.acid_exposure_fraction * 0.7) + (1.0 - self.mucosal_integrity) * 0.8; + let target_drive = (0.2 + dryness + irritation).clamp(0.05, 1.0); + self.swallow_drive = Self::approach(self.swallow_drive, target_drive, 0.35, dt_seconds); + let interval_target = (35.0 - 18.0 * self.swallow_drive).clamp(4.5, 55.0); + self.swallow_interval_target_s = Self::approach( + self.swallow_interval_target_s, + interval_target, + 0.5, + dt_seconds, + ); + let vagal_target = + (0.6 + 0.25 * self.swallow_drive - 0.15 * self.acid_exposure_fraction).clamp(0.35, 0.9); + self.vagal_tone = Self::approach(self.vagal_tone, vagal_target, 0.3, dt_seconds); + let strength_target = (0.85 + 0.5 * (self.vagal_tone - 0.6)).clamp(0.5, 1.3); + self.peristaltic_strength = + Self::approach(self.peristaltic_strength, strength_target, 0.4, dt_seconds); + } + + fn update_sphincter_tones(&mut self, dt_seconds: f32) { + let base_les = + (0.7 + 0.18 * self.vagal_tone - 0.25 * self.acid_exposure_fraction).clamp(0.25, 0.98); + let base_ues = (0.82 + 0.1 * self.vagal_tone - 0.1 * self.swallow_drive).clamp(0.4, 0.98); + let (les_modifier, ues_modifier) = match self.stage { + EsophagealStage::Idle => (0.0, 0.0), + EsophagealStage::SwallowInitiation => (-0.4, -0.5), + EsophagealStage::PrimaryPeristalsis => (-0.25, -0.4), + EsophagealStage::SecondaryPeristalsis => (-0.18, -0.2), + EsophagealStage::Clearing => (-0.1, -0.1), + EsophagealStage::RefluxExposure => (-0.35, 0.0), + }; + let les_target = (base_les + les_modifier).clamp(0.1, 0.98); + let ues_target = (base_ues + ues_modifier).clamp(0.05, 0.98); + self.lower_sphincter_tone = + Self::approach(self.lower_sphincter_tone, les_target, 1.4, dt_seconds).clamp(0.05, 1.0); + self.upper_sphincter_tone = + Self::approach(self.upper_sphincter_tone, ues_target, 2.0, dt_seconds).clamp(0.05, 1.0); + } + + fn update_hiatal_gradient(&mut self, dt_seconds: f32) { + let target = (BASE_HIATAL_PRESSURE_CM_H2O + + 2.0 * (self.swallow_drive - 0.3) + + if matches!(self.stage, EsophagealStage::RefluxExposure) { + 1.0 + } else { + 0.0 + }) + .clamp(3.0, 18.0); + self.hiatal_pressure_gradient_cm_h2o = Self::approach( + self.hiatal_pressure_gradient_cm_h2o, + target, + 0.25, + dt_seconds, + ); + } + + fn handle_stage(&mut self, dt_seconds: f32) { + match self.stage { + EsophagealStage::Idle => { + self.peristaltic_progress_cm = + Self::approach(self.peristaltic_progress_cm, 0.0, 10.0, dt_seconds); + self.bolus_volume_ml = Self::approach(self.bolus_volume_ml, 0.0, 6.0, dt_seconds); + if self.acid_exposure_fraction > 0.35 { + self.transition_stage(EsophagealStage::RefluxExposure); + } + } + EsophagealStage::SwallowInitiation => { + if self.time_in_stage_s > 0.28 { + self.transition_stage(EsophagealStage::PrimaryPeristalsis); + } + } + EsophagealStage::PrimaryPeristalsis => { + let speed = PRIMARY_WAVE_SPEED_CM_S * self.peristaltic_strength.clamp(0.4, 1.5); + self.peristaltic_progress_cm += speed * dt_seconds; + self.bolus_volume_ml = (self.bolus_volume_ml - dt_seconds * (speed * 0.8)).max(0.6); + if self.peristaltic_progress_cm >= ESOPHAGEAL_LENGTH_CM { + if self.bolus_volume_ml > 1.2 { + self.transition_stage(EsophagealStage::SecondaryPeristalsis); + } else { + self.transition_stage(EsophagealStage::Clearing); + } + } + } + EsophagealStage::SecondaryPeristalsis => { + let speed = SECONDARY_WAVE_SPEED_CM_S * self.peristaltic_strength.clamp(0.4, 1.4); + self.peristaltic_progress_cm += speed * dt_seconds; + self.bolus_volume_ml = (self.bolus_volume_ml - dt_seconds * (speed * 0.7)).max(0.3); + if self.time_in_stage_s > 6.0 || self.bolus_volume_ml <= 0.4 { + self.transition_stage(EsophagealStage::Clearing); + } + } + EsophagealStage::Clearing => { + self.bolus_volume_ml = Self::approach(self.bolus_volume_ml, 0.0, 4.0, dt_seconds); + if self.time_in_stage_s > 2.0 { + self.transition_stage(EsophagealStage::Idle); + } + } + EsophagealStage::RefluxExposure => { + self.peristaltic_progress_cm = + Self::approach(self.peristaltic_progress_cm, 0.0, 5.0, dt_seconds); + if self.acid_exposure_fraction < 0.12 { + self.transition_stage(EsophagealStage::Idle); + } + } + } + } + + fn update_acid_balance(&mut self, dt_seconds: f32) { + let reflux_propensity = (self.hiatal_pressure_gradient_cm_h2o - 4.0).max(0.0) + * (1.0 - self.lower_sphincter_tone); + let reflux_influx = reflux_propensity * 0.012 * dt_seconds; + if reflux_influx > 0.0 { + self.acid_exposure_fraction = + (self.acid_exposure_fraction + reflux_influx).clamp(0.0, 1.2); + if self.acid_exposure_fraction > 0.25 { + self.transition_stage(EsophagealStage::RefluxExposure); + } + } + let saliva_clearance = + (self.saliva_buffer_ml * 0.015 + self.peristaltic_strength * 0.01) * dt_seconds; + self.acid_exposure_fraction = (self.acid_exposure_fraction - saliva_clearance).max(0.0); + let mucosal_damage = (self.acid_exposure_fraction - 0.18).max(0.0) * 0.006 * dt_seconds; + let mucosal_healing = (self.saliva_buffer_ml * 0.003 + 0.0025) * dt_seconds; + self.mucosal_integrity = + (self.mucosal_integrity - mucosal_damage + mucosal_healing).clamp(0.55, 1.02); + let target_reflux_rate = (reflux_propensity * 11.0).clamp(0.0, 30.0); + self.reflux_events_per_hour = Self::approach( + self.reflux_events_per_hour, + target_reflux_rate, + 0.12, + dt_seconds, + ); + let acid_drop = (self.acid_exposure_fraction * 4.5).clamp(0.0, 6.5); + let base_recovery = (self.saliva_buffer_ml * 0.18).clamp(0.0, 3.0); + self.luminal_ph = (6.5 - acid_drop + base_recovery).clamp(1.0, 7.2); + self.saliva_buffer_ml = (self.saliva_buffer_ml - dt_seconds * 1.0).max(0.0); + } } impl Organ for Esophagus { @@ -24,11 +257,39 @@ impl Organ for Esophagus { fn organ_type(&self) -> OrganType { self.info.kind() } - fn update(&mut self, _dt_seconds: f32) { - self.reflux = self.reflux.min(100); + fn update(&mut self, dt_seconds: f32) { + if dt_seconds <= 0.0 { + return; + } + + self.time_since_last_swallow_s += dt_seconds; + self.time_in_stage_s += dt_seconds; + + self.update_swallow_drive(dt_seconds); + self.update_hiatal_gradient(dt_seconds); + self.update_sphincter_tones(dt_seconds); + + if self.time_since_last_swallow_s >= self.swallow_interval_target_s + && matches!( + self.stage, + EsophagealStage::Idle | EsophagealStage::RefluxExposure + ) + { + self.start_primary_swallow(); + } + + self.handle_stage(dt_seconds); + self.update_acid_balance(dt_seconds); } fn summary(&self) -> String { - format!("Esophagus[id={}, reflux={}]", self.id(), self.reflux) + format!( + "Esophagus[id={}, stage={:?}, pH={:.1}, LES={:.0}%, reflux≈{:.1}/h]", + self.id(), + self.stage, + self.luminal_ph, + self.lower_sphincter_tone * 100.0, + self.reflux_events_per_hour + ) } fn as_any(&self) -> &dyn core::any::Any { self diff --git a/src/organs/gallbladder.rs b/src/organs/gallbladder.rs index 7e1bc7e..94196f7 100644 --- a/src/organs/gallbladder.rs +++ b/src/organs/gallbladder.rs @@ -1,20 +1,224 @@ use super::{Organ, OrganInfo}; use crate::types::OrganType; +const DEFAULT_CAPACITY_ML: f32 = 50.0; +const RESIDUAL_VOLUME_ML: f32 = 8.0; + +/// Functional state of the gallbladder during the interdigestive and post-prandial cycle. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum GallbladderPhase { + Filling, + Primed, + Contraction, + Expulsion, + Recovery, +} + #[derive(Debug, Clone)] pub struct Gallbladder { info: OrganInfo, - /// Bile volume ml - pub bile_ml: f32, + /// Stored bile volume (ml). + pub bile_volume_ml: f32, + /// Bile acid concentration (mmol/L). + pub bile_acid_concentration_mmol_l: f32, + /// Cholesterol saturation index (dimensionless; >1 predisposes to stones). + pub cholesterol_saturation_index: f32, + /// Flow of bile exiting via the cystic duct/common bile duct (ml/min). + pub bile_outflow_ml_per_min: f32, + /// Tone of the sphincter of Oddi (0..=1, higher = more closed). + pub sphincter_of_oddi_tone: f32, + /// Circulating cholecystokinin level (ng/mL proxy). + pub cck_level_ng_ml: f32, + /// Hepatic bile inflow (ml/min) delivered to the gallbladder when filling. + pub hepatic_bile_flow_ml_per_min: f32, + /// Total bile-acid pool currently stored in the gallbladder (mmol). + pub bile_acid_pool_mmol: f32, + /// Efficiency of enterohepatic recycling (0..=1). + pub bile_salt_recycling_efficiency: f32, + /// Vagal tone supporting coordinated contraction (0..=1). + pub vagal_tone: f32, + /// Fractional mucosal absorption rate. + pub mucosal_absorption_fraction: f32, + /// Current state in the contraction cycle. + pub phase: GallbladderPhase, + /// External or simulated meal stimulus (0..=1). + pub meal_signal: f32, + /// Internal fasting-driven feed-forward signal (0..=1). + pub internal_meal_drive: f32, + time_in_phase_s: f32, + fasting_clock_s: f32, + target_meal_interval_s: f32, + /// Stone-forming propensity index (0..≈2). + pub gallstone_nucleation_index: f32, } impl Gallbladder { pub fn new(id: impl Into) -> Self { Self { info: OrganInfo::new(id, OrganType::Gallbladder), - bile_ml: 30.0, + bile_volume_ml: 35.0, + bile_acid_concentration_mmol_l: 65.0, + cholesterol_saturation_index: 0.9, + bile_outflow_ml_per_min: 0.05, + sphincter_of_oddi_tone: 0.85, + cck_level_ng_ml: 0.2, + hepatic_bile_flow_ml_per_min: 0.55, + bile_acid_pool_mmol: 2.6, + bile_salt_recycling_efficiency: 0.92, + vagal_tone: 0.58, + mucosal_absorption_fraction: 0.03, + phase: GallbladderPhase::Filling, + meal_signal: 0.1, + internal_meal_drive: 0.12, + time_in_phase_s: 0.0, + fasting_clock_s: 0.0, + target_meal_interval_s: 4.8 * 3600.0, + gallstone_nucleation_index: 0.2, } } + + fn approach(current: f32, target: f32, rate_per_second: f32, dt_seconds: f32) -> f32 { + let rate = rate_per_second.max(0.0); + if rate == 0.0 || dt_seconds <= 0.0 { + return current; + } + let delta = target - current; + let max_step = rate * dt_seconds; + if delta > max_step { + current + max_step + } else if delta < -max_step { + current - max_step + } else { + target + } + } + + fn transition_phase(&mut self, phase: GallbladderPhase) { + if self.phase != phase { + self.phase = phase; + self.time_in_phase_s = 0.0; + } + } + + fn update_meal_drives(&mut self, dt_seconds: f32) { + self.fasting_clock_s += dt_seconds; + if self.fasting_clock_s >= self.target_meal_interval_s { + self.internal_meal_drive = 1.0; + self.fasting_clock_s = 0.0; + self.target_meal_interval_s = (4.0 + 1.5 * self.vagal_tone) * 3600.0; + } + self.internal_meal_drive = (self.internal_meal_drive - dt_seconds / 900.0).clamp(0.0, 1.0); + // Allow external stimuli to decay gently if not continuously refreshed. + self.meal_signal = Self::approach(self.meal_signal, 0.1, 0.05, dt_seconds); + let stimulus = self.meal_signal.max(self.internal_meal_drive); + let cck_target = (0.15 + 1.6 * stimulus).clamp(0.1, 2.5); + self.cck_level_ng_ml = + Self::approach(self.cck_level_ng_ml, cck_target, 0.8, dt_seconds).clamp(0.05, 3.0); + let vagal_target = (0.55 + 0.35 * stimulus).clamp(0.4, 0.92); + self.vagal_tone = Self::approach(self.vagal_tone, vagal_target, 0.35, dt_seconds); + } + + fn update_sphincter_tone(&mut self, dt_seconds: f32) { + let tone_target = match self.phase { + GallbladderPhase::Filling => 0.88, + GallbladderPhase::Primed => 0.75, + GallbladderPhase::Contraction => 0.45, + GallbladderPhase::Expulsion => 0.35, + GallbladderPhase::Recovery => 0.7, + } - 0.18 * (self.cck_level_ng_ml - 0.3).max(0.0); + self.sphincter_of_oddi_tone = Self::approach( + self.sphincter_of_oddi_tone, + tone_target.clamp(0.2, 0.95), + 0.9, + dt_seconds, + ); + } + + fn hepatic_inflow(&self) -> f32 { + match self.phase { + GallbladderPhase::Contraction | GallbladderPhase::Expulsion => { + self.hepatic_bile_flow_ml_per_min * 0.4 + } + _ => self.hepatic_bile_flow_ml_per_min, + } + } + + fn update_bile_pool(&mut self, dt_seconds: f32, outflow_ml: f32) { + let inflow_ml = self.hepatic_inflow() * dt_seconds / 60.0; + let inflow_bile_acids = inflow_ml * 0.075 * self.bile_salt_recycling_efficiency; + let absorption_loss = + self.bile_acid_pool_mmol * self.mucosal_absorption_fraction * dt_seconds / 3600.0; + let pool_after = (self.bile_acid_pool_mmol + inflow_bile_acids - absorption_loss).max(0.5); + let volume_after = (self.bile_volume_ml + inflow_ml - outflow_ml) + .clamp(RESIDUAL_VOLUME_ML, DEFAULT_CAPACITY_ML); + let volume_ratio = if self.bile_volume_ml > 0.0 { + outflow_ml / self.bile_volume_ml + } else { + 0.0 + } + .clamp(0.0, 0.95); + let pool_after = (pool_after * (1.0 - volume_ratio)).max(0.5); + self.bile_volume_ml = volume_after; + self.bile_acid_pool_mmol = pool_after; + let volume_l = (self.bile_volume_ml / 1000.0).max(0.005); + self.bile_acid_concentration_mmol_l = + (self.bile_acid_pool_mmol / volume_l).clamp(20.0, 140.0); + let saturation = (1.0 + 0.32 * (self.bile_volume_ml / DEFAULT_CAPACITY_ML - 0.5) + - 0.45 * (self.bile_acid_concentration_mmol_l / 60.0 - 1.0)) + .clamp(0.6, 1.6); + self.cholesterol_saturation_index = saturation; + } + + fn handle_phase(&mut self, _dt_seconds: f32) { + match self.phase { + GallbladderPhase::Filling => { + self.bile_outflow_ml_per_min = 0.05 * (1.0 - self.sphincter_of_oddi_tone); + if (self.cck_level_ng_ml > 0.35 && self.bile_volume_ml > DEFAULT_CAPACITY_ML * 0.55) + || self.mucosal_absorption_fraction < 0.02 + { + self.transition_phase(GallbladderPhase::Primed); + } + } + GallbladderPhase::Primed => { + self.bile_outflow_ml_per_min = 0.2 + 1.5 * (self.cck_level_ng_ml - 0.3).max(0.0); + if self.time_in_phase_s > 60.0 || self.cck_level_ng_ml > 0.6 { + self.transition_phase(GallbladderPhase::Contraction); + } + } + GallbladderPhase::Contraction => { + self.bile_outflow_ml_per_min = (2.5 + + 8.0 * (self.cck_level_ng_ml - 0.3).max(0.0) + + 4.5 * (self.vagal_tone - 0.5).max(0.0)) + * (1.0 - self.sphincter_of_oddi_tone); + if self.bile_volume_ml < DEFAULT_CAPACITY_ML * 0.3 { + self.transition_phase(GallbladderPhase::Expulsion); + } + } + GallbladderPhase::Expulsion => { + self.bile_outflow_ml_per_min = + (1.2 + 6.0 * (1.0 - self.sphincter_of_oddi_tone)).clamp(0.6, 12.0); + if self.bile_volume_ml <= RESIDUAL_VOLUME_ML + 1.0 || self.time_in_phase_s > 180.0 { + self.transition_phase(GallbladderPhase::Recovery); + } + } + GallbladderPhase::Recovery => { + self.bile_outflow_ml_per_min = 0.1 * (1.0 - self.sphincter_of_oddi_tone); + if self.cck_level_ng_ml < 0.25 && self.time_in_phase_s > 120.0 { + self.transition_phase(GallbladderPhase::Filling); + } + } + } + } + + fn update_gallstone_index(&mut self, dt_seconds: f32) { + let stasis = (1.0 - (self.bile_outflow_ml_per_min / 8.0).clamp(0.0, 1.0)).max(0.0); + let supersaturation = (self.cholesterol_saturation_index - 1.0).max(0.0); + let volume_factor = (self.bile_volume_ml / DEFAULT_CAPACITY_ML - 0.6).max(0.0); + let target = + (0.3 + 1.8 * supersaturation * (0.6 + stasis) + 0.5 * volume_factor).clamp(0.0, 2.2); + self.gallstone_nucleation_index = + Self::approach(self.gallstone_nucleation_index, target, 0.15, dt_seconds); + } } impl Organ for Gallbladder { @@ -24,9 +228,48 @@ impl Organ for Gallbladder { fn organ_type(&self) -> OrganType { self.info.kind() } - fn update(&mut self, _dt_seconds: f32) {} + fn update(&mut self, dt_seconds: f32) { + if dt_seconds <= 0.0 { + return; + } + + self.time_in_phase_s += dt_seconds; + + self.update_meal_drives(dt_seconds); + self.update_sphincter_tone(dt_seconds); + self.handle_phase(dt_seconds); + + let outflow_ml = (self.bile_outflow_ml_per_min * dt_seconds / 60.0).clamp(0.0, 25.0); + self.update_bile_pool(dt_seconds, outflow_ml); + + // Adjust mucosal absorption with concentration and phase. + let absorption_target = match self.phase { + GallbladderPhase::Filling => 0.035, + GallbladderPhase::Primed => 0.03, + GallbladderPhase::Contraction => 0.02, + GallbladderPhase::Expulsion => 0.015, + GallbladderPhase::Recovery => 0.028, + } * (self.bile_acid_concentration_mmol_l / 60.0).clamp(0.7, 1.4); + self.mucosal_absorption_fraction = Self::approach( + self.mucosal_absorption_fraction, + absorption_target, + 0.2, + dt_seconds, + ) + .clamp(0.01, 0.05); + + self.update_gallstone_index(dt_seconds); + } fn summary(&self) -> String { - format!("Gallbladder[id={}, bile={:.0} ml]", self.id(), self.bile_ml) + format!( + "Gallbladder[id={}, phase={:?}, vol={:.0}/{:.0} ml, bileAcid={:.0} mmol/L, CSI={:.2}]", + self.id(), + self.phase, + self.bile_volume_ml, + DEFAULT_CAPACITY_ML, + self.bile_acid_concentration_mmol_l, + self.cholesterol_saturation_index + ) } fn as_any(&self) -> &dyn core::any::Any { self diff --git a/src/organs/heart.rs b/src/organs/heart.rs index 6e94c24..1489ae9 100644 --- a/src/organs/heart.rs +++ b/src/organs/heart.rs @@ -1,18 +1,71 @@ use super::{Organ, OrganInfo}; use crate::types::{BloodPressure, OrganType}; -/// Cardiac model with simple rate and arterial pressure coupling. +const BASE_SV_ML: f32 = 70.0; +const BASE_SVR_MMHG_MIN_PER_L: f32 = 17.0; +const BAROREFLEX_SET_POINT_MMHG: f32 = 93.0; + +/// Rhythm archetypes representing dominant autonomic/conduction control of the heart. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum CardiacRhythmState { + Sinus, + SympatheticDrive, + ParasympatheticDominance, + CompensatoryTachycardia, + Arrhythmic, +} + +/// Cardiac pump model featuring autonomic reflexes and hemodynamic coupling. #[derive(Debug, Clone)] pub struct Heart { info: OrganInfo, - /// Heart rate in beats per minute. + /// Heart rate (beats per minute). pub heart_rate_bpm: f32, /// Arterial blood pressure snapshot. pub arterial_bp: BloodPressure, - /// ECG lead count configured for this heart. + /// Number of ECG leads configured for monitoring. pub leads: u8, - /// Simplified arrhythmia flag; increases HR variability. + /// Allow external systems to force arrhythmic behavior. pub arrhythmia: bool, + /// Stroke volume (ml/beat). + pub stroke_volume_ml: f32, + /// Cardiac output (L/min). + pub cardiac_output_l_min: f32, + /// End-diastolic volume (ml). + pub end_diastolic_volume_ml: f32, + /// End-systolic volume (ml). + pub end_systolic_volume_ml: f32, + /// Ejection fraction (0..=1). + pub ejection_fraction: f32, + /// Contractility index (dimensionless relative to baseline 1.0). + pub contractility_index: f32, + /// Preload expressed as estimated left-ventricular end-diastolic pressure (mmHg). + pub preload_mm_hg: f32, + /// Afterload (mmHg) approximated from systemic vascular resistance. + pub afterload_mm_hg: f32, + /// Systemic vascular resistance (mmHg·min/L). + pub systemic_vascular_resistance: f32, + /// Venous return (L/min). + pub venous_return_l_min: f32, + /// Sinoatrial node intrinsic rate (bpm). + pub sa_node_rate_bpm: f32, + /// Atrioventricular conduction delay (ms). + pub av_delay_ms: f32, + /// Autonomic tone (-1 parasympathetic, +1 sympathetic). + pub autonomic_tone: f32, + /// Current rhythm classification. + pub rhythm_state: CardiacRhythmState, + /// Arrhythmia burden (0..=1). + pub arrhythmia_burden: f32, + /// Myocardial oxygen demand (mL O2/beat scaled). + pub myocardial_oxygen_demand: f32, + /// Myocardial oxygen supply proxy (mL O2/beat scaled). + pub myocardial_oxygen_supply: f32, + /// Coronary perfusion pressure (mmHg). + pub coronary_perfusion_mm_hg: f32, + /// Stroke work (J per beat). + pub stroke_work_joule: f32, + time_in_state_s: f32, } impl Heart { @@ -20,12 +73,189 @@ impl Heart { pub fn new(id: impl Into, leads: u8) -> Self { Self { info: OrganInfo::new(id, OrganType::Heart), - heart_rate_bpm: 70.0, + heart_rate_bpm: 72.0, arterial_bp: BloodPressure::default(), leads, arrhythmia: false, + stroke_volume_ml: BASE_SV_ML, + cardiac_output_l_min: 5.0, + end_diastolic_volume_ml: 120.0, + end_systolic_volume_ml: 50.0, + ejection_fraction: 0.58, + contractility_index: 1.0, + preload_mm_hg: 8.0, + afterload_mm_hg: 85.0, + systemic_vascular_resistance: BASE_SVR_MMHG_MIN_PER_L, + venous_return_l_min: 5.0, + sa_node_rate_bpm: 72.0, + av_delay_ms: 160.0, + autonomic_tone: 0.0, + rhythm_state: CardiacRhythmState::Sinus, + arrhythmia_burden: 0.05, + myocardial_oxygen_demand: 9.0, + myocardial_oxygen_supply: 9.5, + coronary_perfusion_mm_hg: 70.0, + stroke_work_joule: 1.1, + time_in_state_s: 0.0, } } + + fn approach(current: f32, target: f32, rate_per_second: f32, dt_seconds: f32) -> f32 { + let rate = rate_per_second.max(0.0); + if rate == 0.0 || dt_seconds <= 0.0 { + return current; + } + let delta = target - current; + let max_step = rate * dt_seconds; + if delta > max_step { + current + max_step + } else if delta < -max_step { + current - max_step + } else { + target + } + } + + fn mean_arterial_pressure(&self) -> f32 { + let sys = self.arterial_bp.systolic as f32; + let dia = self.arterial_bp.diastolic as f32; + dia + (sys - dia) / 3.0 + } + + fn update_autonomic_state(&mut self, dt_seconds: f32) { + let map_error = self.mean_arterial_pressure() - BAROREFLEX_SET_POINT_MMHG; + let autonomic_target = (-map_error / 30.0).clamp(-1.0, 1.0); + self.autonomic_tone = + Self::approach(self.autonomic_tone, autonomic_target, 0.8, dt_seconds); + self.sa_node_rate_bpm = Self::approach( + self.sa_node_rate_bpm, + (70.0 + 45.0 * self.autonomic_tone).clamp(45.0, 150.0), + 1.5, + dt_seconds, + ); + self.av_delay_ms = Self::approach( + self.av_delay_ms, + (170.0 - 40.0 * self.autonomic_tone).clamp(110.0, 240.0), + 5.0, + dt_seconds, + ); + let svr_target = (BASE_SVR_MMHG_MIN_PER_L - 5.5 * self.autonomic_tone).clamp(10.0, 26.0); + self.systemic_vascular_resistance = Self::approach( + self.systemic_vascular_resistance, + svr_target, + 0.3, + dt_seconds, + ); + } + + fn determine_rhythm_state(&mut self) { + let arrhythmic = self.arrhythmia || self.arrhythmia_burden > 0.45; + self.rhythm_state = if arrhythmic { + CardiacRhythmState::Arrhythmic + } else if self.autonomic_tone > 0.5 { + CardiacRhythmState::SympatheticDrive + } else if self.autonomic_tone < -0.4 { + CardiacRhythmState::ParasympatheticDominance + } else if self.venous_return_l_min < 4.2 && self.mean_arterial_pressure() < 75.0 { + CardiacRhythmState::CompensatoryTachycardia + } else { + CardiacRhythmState::Sinus + }; + } + + fn update_rate_and_contractility(&mut self, dt_seconds: f32) { + let mut rate_target = self.sa_node_rate_bpm; + match self.rhythm_state { + CardiacRhythmState::SympatheticDrive => rate_target += 18.0, + CardiacRhythmState::ParasympatheticDominance => rate_target -= 12.0, + CardiacRhythmState::CompensatoryTachycardia => rate_target += 22.0, + CardiacRhythmState::Arrhythmic => rate_target += 8.0, + CardiacRhythmState::Sinus => {} + } + rate_target += 8.0 * (self.arrhythmia_burden - 0.2).max(0.0); + self.heart_rate_bpm = Self::approach( + self.heart_rate_bpm, + rate_target.clamp(38.0, 190.0), + 1.2, + dt_seconds, + ); + let demand_scale = (self.heart_rate_bpm / 70.0).clamp(0.6, 2.0); + let afterload_penalty = (self.afterload_mm_hg - 90.0).max(0.0) / 120.0; + let contractility_target = (1.05 + 0.35 * self.autonomic_tone + - afterload_penalty + - (self.arrhythmia_burden * 0.2)) + .clamp(0.5, 1.6); + self.contractility_index = Self::approach( + self.contractility_index, + contractility_target, + 0.8, + dt_seconds, + ); + self.myocardial_oxygen_demand = (8.5 * demand_scale + + 0.6 * self.contractility_index * (self.afterload_mm_hg / 80.0)) + .clamp(4.0, 20.0); + } + + fn update_volumes_and_output(&mut self, dt_seconds: f32) { + self.preload_mm_hg = Self::approach( + self.preload_mm_hg, + (6.5 + 1.2 * (self.venous_return_l_min - 4.5)).clamp(4.0, 18.0), + 0.6, + dt_seconds, + ); + let edv_target = (90.0 + 6.5 * self.preload_mm_hg).clamp(80.0, 210.0); + self.end_diastolic_volume_ml = + Self::approach(self.end_diastolic_volume_ml, edv_target, 1.1, dt_seconds); + let elastance = 0.22 + 0.25 * self.contractility_index; + let esv_target = (self.end_diastolic_volume_ml * (1.0 - elastance)).clamp(30.0, 120.0); + self.end_systolic_volume_ml = + Self::approach(self.end_systolic_volume_ml, esv_target, 1.6, dt_seconds); + self.stroke_volume_ml = + (self.end_diastolic_volume_ml - self.end_systolic_volume_ml).clamp(25.0, 130.0); + self.ejection_fraction = + (self.stroke_volume_ml / self.end_diastolic_volume_ml.max(1.0)).clamp(0.2, 0.85); + self.cardiac_output_l_min = + (self.stroke_volume_ml * self.heart_rate_bpm / 1000.0).clamp(2.0, 12.0); + self.venous_return_l_min = Self::approach( + self.venous_return_l_min, + self.cardiac_output_l_min, + 0.4, + dt_seconds, + ); + self.afterload_mm_hg = Self::approach( + self.afterload_mm_hg, + (self.systemic_vascular_resistance * self.cardiac_output_l_min).clamp(60.0, 160.0), + 0.5, + dt_seconds, + ); + let map_target = self.cardiac_output_l_min * self.systemic_vascular_resistance; + let pulse_pressure = (self.stroke_volume_ml / BASE_SV_ML).clamp(0.6, 2.0) * 40.0; + let systolic = (map_target + pulse_pressure / 2.0).clamp(80.0, 220.0); + let diastolic = (map_target - pulse_pressure / 2.5).clamp(40.0, systolic - 5.0); + self.arterial_bp.systolic = systolic.round() as u16; + self.arterial_bp.diastolic = diastolic.round() as u16; + self.coronary_perfusion_mm_hg = + (self.arterial_bp.diastolic as f32 - self.preload_mm_hg).clamp(20.0, 120.0); + self.myocardial_oxygen_supply = (9.5 + 0.08 * self.coronary_perfusion_mm_hg + - 0.5 * (self.heart_rate_bpm - 70.0) / 40.0) + .clamp(4.0, 20.0); + self.stroke_work_joule = (self.stroke_volume_ml / 1000.0) + * (self.mean_arterial_pressure() - 5.0).max(0.0) + * 0.133; + } + + fn update_arrhythmia_burden(&mut self, dt_seconds: f32) { + let supply_demand_ratio = + (self.myocardial_oxygen_supply / self.myocardial_oxygen_demand).clamp(0.4, 1.6); + let mismatch = (1.0 - supply_demand_ratio).max(0.0); + let target = if self.arrhythmia { + 0.7 + } else { + 0.2 + 0.6 * mismatch + 0.2 * (self.autonomic_tone - 0.5).max(0.0) + } + .clamp(0.05, 0.95); + self.arrhythmia_burden = Self::approach(self.arrhythmia_burden, target, 0.3, dt_seconds); + } } impl Organ for Heart { @@ -36,26 +266,27 @@ impl Organ for Heart { self.info.kind() } fn update(&mut self, dt_seconds: f32) { - let dt = dt_seconds.clamp(0.0, 10.0); - let target = 70.0f32; - let mut diff = target - self.heart_rate_bpm; - if self.arrhythmia { - // add variability - diff += self.heart_rate_bpm.sin() * 5.0; + if dt_seconds <= 0.0 { + return; } - self.heart_rate_bpm += 0.1 * diff * (dt / 1.0); - // crude BP coupling to HR - let sys = (90.0 + 0.5 * self.heart_rate_bpm).clamp(80.0, 220.0) as u16; - let dia = (60.0 + 0.3 * self.heart_rate_bpm).clamp(40.0, 140.0) as u16; - self.arterial_bp.systolic = sys; - self.arterial_bp.diastolic = dia.min(sys.saturating_sub(10)); + + self.time_in_state_s += dt_seconds; + + self.update_autonomic_state(dt_seconds); + self.determine_rhythm_state(); + self.update_rate_and_contractility(dt_seconds); + self.update_volumes_and_output(dt_seconds); + self.update_arrhythmia_burden(dt_seconds); } fn summary(&self) -> String { format!( - "Heart[id={}, leads={}, HR={:.1} bpm, BP={}/{} mmHg]", + "Heart[id={}, leads={}, rhythm={:?}, HR={:.0} bpm, CO={:.1} L/min, EF={:.0}%, BP={}/{}]", self.id(), self.leads, + self.rhythm_state, self.heart_rate_bpm, + self.cardiac_output_l_min, + self.ejection_fraction * 100.0, self.arterial_bp.systolic, self.arterial_bp.diastolic ) diff --git a/src/organs/intestines.rs b/src/organs/intestines.rs index 76676f3..5703b74 100644 --- a/src/organs/intestines.rs +++ b/src/organs/intestines.rs @@ -1,22 +1,293 @@ use super::{Organ, OrganInfo}; use crate::types::OrganType; +/// Predominant motility/functional state of the small and large intestine. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum IntestinalPhase { + Interdigestive, + FedProcessing, + MigratingMotorComplex, + IlealBrake, + Dysmotility, +} + #[derive(Debug, Clone)] pub struct Intestines { info: OrganInfo, - /// Nutrient absorption rate 0..=100 - pub absorption: u8, - pub peristalsis: bool, + /// Carbohydrate absorption (g/hour). + pub carbohydrate_absorption_g_per_h: f32, + /// Fat absorption (g/hour). + pub fat_absorption_g_per_h: f32, + /// Protein absorption (g/hour). + pub protein_absorption_g_per_h: f32, + /// Electrolyte reclamation (mmol/min). + pub electrolyte_absorption_mmol_min: f32, + /// Water reabsorption rate (ml/min). + pub water_reabsorption_ml_min: f32, + /// Integrated motility index (0..=1) dominated by peristaltic waves. + pub motility_index: f32, + /// Segmentation/segmental contractions index (0..=1). + pub segmentation_index: f32, + /// Migrating motor complex activity (0..=1). + pub mmc_activity: f32, + /// Current functional phase. + pub phase: IntestinalPhase, + /// Luminal volume of chyme (ml). + pub lumen_volume_ml: f32, + /// Luminal pH. + pub chyme_ph: f32, + /// Fraction of bile acids reclaimed in the terminal ileum (0..=1). + pub bile_acid_recirculation_fraction: f32, + /// Microbiome balance (0..=1, >0.5 reflects eubiosis). + pub microbiome_balance: f32, + /// Short-chain fatty acids produced (mmol). + pub short_chain_fatty_acids_mmol: f32, + /// Mucosal integrity (0..=1). + pub mucosal_integrity: f32, + /// Inflammatory index (0..=1). + pub inflammation_index: f32, + /// Motilin level (pg/mL proxy) driving MMC. + pub hormone_motilin: f32, + /// GLP-1 level influencing ileal brake (pmol/L proxy). + pub hormone_glp1: f32, + /// Enteric nervous system tone (0..=1). + pub enteric_tone: f32, + /// Pending nutrient energy load within the lumen (kcal). + pub nutrient_energy_kcal: f32, + /// Fermentable fiber load (g). + pub fiber_load_g: f32, + time_in_phase_s: f32, + feeding_clock_s: f32, + target_feed_interval_s: f32, } impl Intestines { pub fn new(id: impl Into) -> Self { Self { info: OrganInfo::new(id, OrganType::Intestines), - absorption: 80, - peristalsis: true, + carbohydrate_absorption_g_per_h: 45.0, + fat_absorption_g_per_h: 12.0, + protein_absorption_g_per_h: 18.0, + electrolyte_absorption_mmol_min: 2.5, + water_reabsorption_ml_min: 12.0, + motility_index: 0.55, + segmentation_index: 0.45, + mmc_activity: 0.3, + phase: IntestinalPhase::Interdigestive, + lumen_volume_ml: 350.0, + chyme_ph: 6.6, + bile_acid_recirculation_fraction: 0.92, + microbiome_balance: 0.62, + short_chain_fatty_acids_mmol: 22.0, + mucosal_integrity: 0.93, + inflammation_index: 0.12, + hormone_motilin: 110.0, + hormone_glp1: 8.0, + enteric_tone: 0.55, + nutrient_energy_kcal: 40.0, + fiber_load_g: 6.0, + time_in_phase_s: 0.0, + feeding_clock_s: 0.0, + target_feed_interval_s: 3.8 * 3600.0, } } + + fn approach(current: f32, target: f32, rate_per_second: f32, dt_seconds: f32) -> f32 { + let rate = rate_per_second.max(0.0); + if rate == 0.0 || dt_seconds <= 0.0 { + return current; + } + let delta = target - current; + let max_step = rate * dt_seconds; + if delta > max_step { + current + max_step + } else if delta < -max_step { + current - max_step + } else { + target + } + } + + fn update_internal_feeding(&mut self, dt_seconds: f32) { + self.feeding_clock_s += dt_seconds; + if self.feeding_clock_s >= self.target_feed_interval_s { + self.nutrient_energy_kcal += 410.0; + self.fiber_load_g += 4.0; + self.lumen_volume_ml = (self.lumen_volume_ml + 200.0).clamp(150.0, 800.0); + self.phase = IntestinalPhase::FedProcessing; + self.time_in_phase_s = 0.0; + self.feeding_clock_s = 0.0; + self.target_feed_interval_s = (3.0 + 1.5 * (1.0 - self.microbiome_balance)) * 3600.0; + } + } + + fn update_enteric_tone(&mut self, dt_seconds: f32) { + let load_factor = (self.nutrient_energy_kcal / 400.0).clamp(0.0, 2.0); + let irritation = (1.0 - self.mucosal_integrity) * 1.2 + self.inflammation_index * 0.8; + let tone_target = (0.5 + 0.35 * load_factor - 0.2 * irritation).clamp(0.2, 0.95); + self.enteric_tone = Self::approach(self.enteric_tone, tone_target, 0.4, dt_seconds); + } + + fn transition_phase(&mut self, phase: IntestinalPhase) { + if self.phase != phase { + self.phase = phase; + self.time_in_phase_s = 0.0; + } + } + + fn update_phase(&mut self) { + match self.phase { + IntestinalPhase::Interdigestive => { + if self.nutrient_energy_kcal > 80.0 { + self.transition_phase(IntestinalPhase::FedProcessing); + } else if self.time_in_phase_s > 90.0 * 60.0 { + self.transition_phase(IntestinalPhase::MigratingMotorComplex); + } + } + IntestinalPhase::FedProcessing => { + if self.nutrient_energy_kcal < 100.0 { + self.transition_phase(IntestinalPhase::IlealBrake); + } + } + IntestinalPhase::MigratingMotorComplex => { + if self.nutrient_energy_kcal > 40.0 { + self.transition_phase(IntestinalPhase::FedProcessing); + } else if self.time_in_phase_s > 120.0 * 60.0 { + self.transition_phase(IntestinalPhase::Interdigestive); + } + } + IntestinalPhase::IlealBrake => { + if self.hormone_glp1 < 6.0 && self.fiber_load_g < 8.0 { + self.transition_phase(IntestinalPhase::Interdigestive); + } else if self.mucosal_integrity < 0.6 { + self.transition_phase(IntestinalPhase::Dysmotility); + } + } + IntestinalPhase::Dysmotility => { + if self.mucosal_integrity > 0.75 && self.inflammation_index < 0.3 { + self.transition_phase(IntestinalPhase::Interdigestive); + } + } + } + } + + fn update_motility(&mut self, dt_seconds: f32) { + let (motility_target, segmentation_target, mmc_target, glp1_target, motilin_target) = + match self.phase { + IntestinalPhase::Interdigestive => (0.45, 0.4, 0.65, 6.0, 150.0), + IntestinalPhase::FedProcessing => (0.72, 0.55, 0.25, 18.0, 80.0), + IntestinalPhase::MigratingMotorComplex => (0.6, 0.35, 0.9, 8.0, 210.0), + IntestinalPhase::IlealBrake => (0.38, 0.6, 0.2, 28.0, 70.0), + IntestinalPhase::Dysmotility => (0.28, 0.35, 0.1, 22.0, 60.0), + }; + self.motility_index = Self::approach(self.motility_index, motility_target, 0.6, dt_seconds); + self.segmentation_index = Self::approach( + self.segmentation_index, + segmentation_target, + 0.5, + dt_seconds, + ); + self.mmc_activity = Self::approach(self.mmc_activity, mmc_target, 0.4, dt_seconds); + self.hormone_glp1 = + Self::approach(self.hormone_glp1, glp1_target, 0.2, dt_seconds).clamp(2.0, 35.0); + self.hormone_motilin = + Self::approach(self.hormone_motilin, motilin_target, 0.8, dt_seconds) + .clamp(40.0, 280.0); + } + + fn update_absorption(&mut self, dt_seconds: f32) { + let motility_effect = + (self.motility_index * 1.1 + self.segmentation_index * 0.6).clamp(0.3, 1.6); + let available = self.nutrient_energy_kcal.max(0.0); + let carbs_available = available * 0.55; + let fat_available = available * 0.28; + let protein_available = available * 0.17; + let carbs_abs_target = (carbs_available / 4.0).clamp(0.0, 90.0) * motility_effect; + let fat_abs_target = (fat_available / 9.0).clamp(0.0, 35.0) * motility_effect; + let protein_abs_target = (protein_available / 4.0).clamp(0.0, 50.0) * motility_effect; + self.carbohydrate_absorption_g_per_h = Self::approach( + self.carbohydrate_absorption_g_per_h, + carbs_abs_target, + 0.3, + dt_seconds, + ); + self.fat_absorption_g_per_h = + Self::approach(self.fat_absorption_g_per_h, fat_abs_target, 0.3, dt_seconds); + self.protein_absorption_g_per_h = Self::approach( + self.protein_absorption_g_per_h, + protein_abs_target, + 0.3, + dt_seconds, + ); + let absorbed_kcal = (self.carbohydrate_absorption_g_per_h * 4.0 + + self.protein_absorption_g_per_h * 4.0 + + self.fat_absorption_g_per_h * 9.0) + * dt_seconds + / 3600.0; + self.nutrient_energy_kcal = (self.nutrient_energy_kcal - absorbed_kcal).max(0.0); + self.electrolyte_absorption_mmol_min = Self::approach( + self.electrolyte_absorption_mmol_min, + (2.0 + 0.8 * self.motility_index + 1.2 * self.segmentation_index).clamp(1.0, 6.0), + 0.2, + dt_seconds, + ); + self.water_reabsorption_ml_min = Self::approach( + self.water_reabsorption_ml_min, + (10.0 + 8.0 * self.electrolyte_absorption_mmol_min / 3.0).clamp(5.0, 30.0), + 0.2, + dt_seconds, + ); + let water_removed = self.water_reabsorption_ml_min * dt_seconds / 60.0; + self.lumen_volume_ml = + (self.lumen_volume_ml + absorbed_kcal * 0.2 - water_removed).clamp(120.0, 800.0); + } + + fn update_microbiome(&mut self, dt_seconds: f32) { + let scfa_generation = + (self.fiber_load_g * 0.12 * self.microbiome_balance) * dt_seconds / 60.0; + self.short_chain_fatty_acids_mmol = + (self.short_chain_fatty_acids_mmol + scfa_generation).clamp(5.0, 80.0); + self.fiber_load_g = (self.fiber_load_g - scfa_generation * 0.45).max(0.0); + let ph_target = (6.6 - 0.015 * self.short_chain_fatty_acids_mmol + + 0.2 * (1.0 - self.bile_acid_recirculation_fraction)) + .clamp(5.8, 7.2); + self.chyme_ph = Self::approach(self.chyme_ph, ph_target, 0.1, dt_seconds); + let microbiome_target = (0.6 + 0.15 * (self.short_chain_fatty_acids_mmol / 30.0) + - 0.25 * self.inflammation_index) + .clamp(0.3, 0.95); + self.microbiome_balance = + Self::approach(self.microbiome_balance, microbiome_target, 0.05, dt_seconds); + } + + fn update_mucosa(&mut self, dt_seconds: f32) { + let irritation = (1.0 - self.chyme_ph / 7.0).max(0.0) + + (self.short_chain_fatty_acids_mmol / 50.0).max(0.0) * 0.2; + let mucosal_target = + (0.95 - 0.25 * irritation + 0.1 * self.microbiome_balance).clamp(0.5, 0.99); + self.mucosal_integrity = + Self::approach(self.mucosal_integrity, mucosal_target, 0.02, dt_seconds); + let inflammation_target = + (0.1 + 0.3 * (1.0 - self.mucosal_integrity) + 0.2 * (1.0 - self.microbiome_balance)) + .clamp(0.05, 0.9); + self.inflammation_index = Self::approach( + self.inflammation_index, + inflammation_target, + 0.03, + dt_seconds, + ); + let bile_target = + (0.9 + 0.15 * self.motility_index - 0.1 * self.glp1_effect()).clamp(0.6, 0.98); + self.bile_acid_recirculation_fraction = Self::approach( + self.bile_acid_recirculation_fraction, + bile_target, + 0.03, + dt_seconds, + ); + } + + fn glp1_effect(&self) -> f32 { + (self.hormone_glp1 / 20.0).clamp(0.0, 1.2) + } } impl Organ for Intestines { @@ -27,18 +298,28 @@ impl Organ for Intestines { self.info.kind() } fn update(&mut self, dt_seconds: f32) { - if self.peristalsis { - // minor oscillation around 80 - let delta = (dt_seconds * 0.1).sin(); - let val = (self.absorption as f32 + delta).clamp(0.0, 100.0); - self.absorption = val as u8; + if dt_seconds <= 0.0 { + return; } + + self.time_in_phase_s += dt_seconds; + + self.update_internal_feeding(dt_seconds); + self.update_enteric_tone(dt_seconds); + self.update_phase(); + self.update_motility(dt_seconds); + self.update_absorption(dt_seconds); + self.update_microbiome(dt_seconds); + self.update_mucosa(dt_seconds); } fn summary(&self) -> String { format!( - "Intestines[id={}, absorption={}]", + "Intestines[id={}, phase={:?}, motility={:.2}, lumen={:.0} ml, pH={:.1}]", self.id(), - self.absorption + self.phase, + self.motility_index, + self.lumen_volume_ml, + self.chyme_ph ) } fn as_any(&self) -> &dyn core::any::Any { diff --git a/src/organs/kidneys.rs b/src/organs/kidneys.rs index b8045ee..f8d592f 100644 --- a/src/organs/kidneys.rs +++ b/src/organs/kidneys.rs @@ -1,23 +1,233 @@ use super::{Organ, OrganInfo}; use crate::types::OrganType; +/// Renal perfusion/autoregulation status. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum RenalAutoregulationState { + Autoregulated, + Hypoperfused, + Hyperperfused, + Obstructed, +} + #[derive(Debug, Clone)] pub struct Kidneys { info: OrganInfo, - /// Filtration rate ml/min + /// Glomerular filtration rate (ml/min). pub gfr: f32, - /// Electrolyte balance index -1..=1 + /// Electrolyte balance index -1..=1 (positive = hypernatremia tendency). pub electrolyte_balance: f32, + /// Renal plasma flow (ml/min). + pub renal_plasma_flow_ml_min: f32, + /// Filtration fraction (GFR/RPF). + pub filtration_fraction: f32, + /// Urine osmolality (mOsm/kg). + pub urine_osmolality_mosm: f32, + /// Urine flow (ml/min). + pub urine_flow_ml_min: f32, + /// Renin release proxy (ng/mL/h relative). + pub renin_release: f32, + /// Aldosterone drive (0..=1). + pub aldosterone_drive: f32, + /// Antidiuretic hormone sensitivity (0..=1). + pub adh_sensitivity: f32, + /// Acid-base compensation index (-1 acidotic .. +1 alkalotic). + pub acid_base_balance: f32, + /// Erythropoietin secretion (IU/day equivalent). + pub erythropoietin_iu_per_day: f32, + /// Sympathetic tone to the juxtaglomerular apparatus (0..=1). + pub sympathetic_tone: f32, + /// Tubular sodium reabsorption fraction (0..=1). + pub tubular_na_reabsorption: f32, + /// Potassium excretion (mmol/min). + pub potassium_excretion_mmol_min: f32, + /// Urea excretion (mg/min). + pub urea_excretion_mg_min: f32, + /// Medullary tonicity (mOsm/kg). + pub medullary_tonicity_mosm: f32, + /// Serum osmolality (mOsm/kg). + pub serum_osmolality_mosm: f32, + /// Estimated plasma volume (L). + pub plasma_volume_l: f32, + /// Current autoregulation state. + pub state: RenalAutoregulationState, + time_in_state_s: f32, } impl Kidneys { pub fn new(id: impl Into) -> Self { Self { info: OrganInfo::new(id, OrganType::Kidneys), - gfr: 100.0, + gfr: 110.0, electrolyte_balance: 0.0, + renal_plasma_flow_ml_min: 600.0, + filtration_fraction: 0.18, + urine_osmolality_mosm: 550.0, + urine_flow_ml_min: 1.2, + renin_release: 0.35, + aldosterone_drive: 0.45, + adh_sensitivity: 0.65, + acid_base_balance: 0.0, + erythropoietin_iu_per_day: 18.0, + sympathetic_tone: 0.35, + tubular_na_reabsorption: 0.99, + potassium_excretion_mmol_min: 0.11, + urea_excretion_mg_min: 550.0, + medullary_tonicity_mosm: 750.0, + serum_osmolality_mosm: 290.0, + plasma_volume_l: 3.1, + state: RenalAutoregulationState::Autoregulated, + time_in_state_s: 0.0, } } + + fn approach(current: f32, target: f32, rate_per_second: f32, dt_seconds: f32) -> f32 { + let rate = rate_per_second.max(0.0); + if rate == 0.0 || dt_seconds <= 0.0 { + return current; + } + let delta = target - current; + let max_step = rate * dt_seconds; + if delta > max_step { + current + max_step + } else if delta < -max_step { + current - max_step + } else { + target + } + } + + fn update_state(&mut self) { + let perfusion_ratio = self.renal_plasma_flow_ml_min / 600.0; + let obstruction = (1.0 - self.urine_flow_ml_min / 1.2).max(0.0) + + (self.medullary_tonicity_mosm - 1000.0).max(0.0) / 800.0; + self.state = if obstruction > 0.6 { + RenalAutoregulationState::Obstructed + } else if perfusion_ratio < 0.75 { + RenalAutoregulationState::Hypoperfused + } else if perfusion_ratio > 1.3 { + RenalAutoregulationState::Hyperperfused + } else { + RenalAutoregulationState::Autoregulated + }; + } + + fn update_perfusion(&mut self, dt_seconds: f32) { + let sympathetic_target = + (0.3 + 0.6 * (1.0 - self.plasma_volume_l / 3.0).max(0.0)).clamp(0.1, 0.95); + self.sympathetic_tone = + Self::approach(self.sympathetic_tone, sympathetic_target, 0.2, dt_seconds); + let rpf_target = (600.0 - 220.0 * self.sympathetic_tone + + 120.0 * (self.serum_osmolality_mosm - 285.0) / 20.0) + .clamp(300.0, 950.0); + self.renal_plasma_flow_ml_min = + Self::approach(self.renal_plasma_flow_ml_min, rpf_target, 3.0, dt_seconds); + let filtration_target = match self.state { + RenalAutoregulationState::Autoregulated => 0.18, + RenalAutoregulationState::Hypoperfused => 0.23, + RenalAutoregulationState::Hyperperfused => 0.16, + RenalAutoregulationState::Obstructed => 0.12, + }; + self.filtration_fraction = + Self::approach(self.filtration_fraction, filtration_target, 0.1, dt_seconds) + .clamp(0.08, 0.3); + self.gfr = (self.renal_plasma_flow_ml_min * self.filtration_fraction).clamp(20.0, 180.0); + } + + fn update_hormonal_axes(&mut self, dt_seconds: f32) { + let pressure_error = (105.0 - self.gfr).max(-40.0); + let osmo_error = (self.serum_osmolality_mosm - 285.0) / 15.0; + let renin_target = + (0.25 + 0.015 * pressure_error + 0.4 * self.sympathetic_tone).clamp(0.05, 1.6); + self.renin_release = Self::approach(self.renin_release, renin_target, 0.3, dt_seconds); + let aldosterone_target = + (0.4 + 0.5 * self.renin_release - 0.3 * self.electrolyte_balance).clamp(0.1, 1.0); + self.aldosterone_drive = + Self::approach(self.aldosterone_drive, aldosterone_target, 0.2, dt_seconds); + let adh_target = + (0.55 + 0.35 * osmo_error + 0.25 * (self.aldosterone_drive - 0.4)).clamp(0.1, 1.1); + self.adh_sensitivity = Self::approach(self.adh_sensitivity, adh_target, 0.25, dt_seconds); + } + + fn update_tubular_handling(&mut self, dt_seconds: f32) { + let sodium_target = match self.state { + RenalAutoregulationState::Hypoperfused => 0.995, + RenalAutoregulationState::Autoregulated => 0.99, + RenalAutoregulationState::Hyperperfused => 0.985, + RenalAutoregulationState::Obstructed => 0.975, + } + 0.01 * (self.aldosterone_drive - 0.5); + self.tubular_na_reabsorption = Self::approach( + self.tubular_na_reabsorption, + sodium_target.clamp(0.94, 0.998), + 0.05, + dt_seconds, + ); + let filtered_nacl = self.gfr * 140.0 / 1000.0; // mmol/min approximated + let na_excreted = filtered_nacl * (1.0 - self.tubular_na_reabsorption); + self.electrolyte_balance = (0.5 - na_excreted / 8.0).clamp(-1.2, 1.2); + self.potassium_excretion_mmol_min = Self::approach( + self.potassium_excretion_mmol_min, + (0.08 + 0.12 * self.aldosterone_drive + 0.04 * self.electrolyte_balance) + .clamp(0.02, 0.3), + 0.1, + dt_seconds, + ); + let osmotic_load = 2.1 * na_excreted + self.potassium_excretion_mmol_min * 1.5; + let adh_effect = (1.3 - self.adh_sensitivity).clamp(0.2, 1.3); + self.urine_flow_ml_min = Self::approach( + self.urine_flow_ml_min, + (osmotic_load / 4.0 * adh_effect).clamp(0.2, 10.0), + 0.4, + dt_seconds, + ); + self.urine_osmolality_mosm = Self::approach( + self.urine_osmolality_mosm, + (550.0 + 220.0 * (self.adh_sensitivity - 0.5) + - 120.0 * (self.urine_flow_ml_min - 1.0) / 4.0) + .clamp(120.0, 1200.0), + 0.6, + dt_seconds, + ); + self.medullary_tonicity_mosm = Self::approach( + self.medullary_tonicity_mosm, + (750.0 + 200.0 * (self.adh_sensitivity - 0.6)).clamp(400.0, 1200.0), + 0.1, + dt_seconds, + ); + } + + fn update_acid_base(&mut self, dt_seconds: f32) { + let acid_target = (0.1 * (self.potassium_excretion_mmol_min - 0.12) + + 0.3 * (self.urine_osmolality_mosm - 600.0) / 400.0) + .clamp(-1.0, 1.0); + self.acid_base_balance = + Self::approach(self.acid_base_balance, -acid_target, 0.1, dt_seconds); + let urea_target = (500.0 + 1.5 * (self.gfr - 110.0)).clamp(200.0, 900.0); + self.urea_excretion_mg_min = + Self::approach(self.urea_excretion_mg_min, urea_target, 0.3, dt_seconds); + self.serum_osmolality_mosm = Self::approach( + self.serum_osmolality_mosm, + (285.0 + 5.0 * self.acid_base_balance - 4.0 * self.urine_flow_ml_min / 5.0) + .clamp(270.0, 310.0), + 0.05, + dt_seconds, + ); + let plasma_target = (3.1 + 0.2 * (self.urine_flow_ml_min - 1.2) + - 0.25 * self.acid_base_balance) + .clamp(2.2, 3.8); + self.plasma_volume_l = + Self::approach(self.plasma_volume_l, plasma_target, 0.04, dt_seconds); + } + + fn update_erythropoietin(&mut self, dt_seconds: f32) { + let epo_target = (18.0 + 40.0 * (0.95 - self.median_oxygenation())).clamp(8.0, 45.0); + self.erythropoietin_iu_per_day = + Self::approach(self.erythropoietin_iu_per_day, epo_target, 0.05, dt_seconds); + } + + fn median_oxygenation(&self) -> f32 { + (self.renal_plasma_flow_ml_min / 600.0).clamp(0.5, 1.3) + } } impl Organ for Kidneys { @@ -27,12 +237,29 @@ impl Organ for Kidneys { fn organ_type(&self) -> OrganType { self.info.kind() } - fn update(&mut self, _dt_seconds: f32) { - self.gfr = self.gfr.clamp(0.0, 200.0); - self.electrolyte_balance = self.electrolyte_balance.clamp(-1.0, 1.0); + fn update(&mut self, dt_seconds: f32) { + if dt_seconds <= 0.0 { + return; + } + + self.time_in_state_s += dt_seconds; + self.update_state(); + self.update_perfusion(dt_seconds); + self.update_state(); + self.update_hormonal_axes(dt_seconds); + self.update_tubular_handling(dt_seconds); + self.update_acid_base(dt_seconds); + self.update_erythropoietin(dt_seconds); } fn summary(&self) -> String { - format!("Kidneys[id={}, GFR={:.0} ml/min]", self.id(), self.gfr) + format!( + "Kidneys[id={}, state={:?}, GFR={:.0} ml/min, urine={:.1} ml/min @ {:.0} mOsm]", + self.id(), + self.state, + self.gfr, + self.urine_flow_ml_min, + self.urine_osmolality_mosm + ) } fn as_any(&self) -> &dyn core::any::Any { self diff --git a/src/organs/liver.rs b/src/organs/liver.rs index 55df4f1..ef15f51 100644 --- a/src/organs/liver.rs +++ b/src/organs/liver.rs @@ -1,15 +1,68 @@ use super::{Organ, OrganInfo}; use crate::types::OrganType; +/// High-level metabolic mode of the liver. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum HepaticState { + Postabsorptive, + FedAnabolic, + FastingCatabolic, + AcutePhaseResponse, + Regenerating, +} + #[derive(Debug, Clone)] pub struct Liver { info: OrganInfo, - /// Detox capacity 0..=100 + /// Detox capacity 0..=100. pub detox: u8, - /// Metabolism index + /// Composite metabolic activity index. pub metabolism: f32, - /// Enzyme production index + /// Microsomal enzyme (CYP) activity index. pub enzymes: f32, + /// Hepatic glycogen store (g). + pub glycogen_store_g: f32, + /// Gluconeogenesis rate (mg/kg/min proxy for 70kg adult). + pub gluconeogenesis_rate: f32, + /// Glycogenolysis rate (g/hour). + pub glycogenolysis_rate_g_per_h: f32, + /// De novo lipogenesis rate (g/hour). + pub lipogenesis_rate_g_per_h: f32, + /// Beta-oxidation rate (g/hour). + pub beta_oxidation_rate_g_per_h: f32, + /// Ammonia clearance (µmol/min). + pub ammonia_clearance_umol_min: f32, + /// Plasma albumin concentration (g/dL). + pub albumin_g_dl: f32, + /// Clotting factor synthesis (% of normal). + pub clotting_factor_synthesis_pct: f32, + /// Bile acid synthesis (mg/min). + pub bile_acid_synthesis_mg_min: f32, + /// Bile secretion (ml/min). + pub bile_secretion_ml_min: f32, + /// Kupffer cell activation (0..=1). + pub kupffer_activation: f32, + /// Acute phase response magnitude (0..=1). + pub acute_phase_response: f32, + /// Hepatic blood flow (L/min). + pub hepatic_blood_flow_l_min: f32, + /// Portal venous pressure (mmHg). + pub portal_pressure_mm_hg: f32, + /// Hepatic fat fraction (%). + pub hepatic_fat_fraction_pct: f32, + /// Insulin signaling intensity (0..=1). + pub insulin_signal: f32, + /// Glucagon signaling intensity (0..=1). + pub glucagon_signal: f32, + /// Cortisol signaling (0..=1). + pub cortisol_signal: f32, + /// Oxidative stress metric (0..=1). + pub oxidative_stress_index: f32, + /// Current metabolic state. + pub state: HepaticState, + time_in_state_s: f32, + feeding_clock_s: f32, + target_meal_interval_s: f32, } impl Liver { @@ -19,8 +72,269 @@ impl Liver { detox: 100, metabolism: 1.0, enzymes: 1.0, + glycogen_store_g: 85.0, + gluconeogenesis_rate: 1.8, + glycogenolysis_rate_g_per_h: 6.0, + lipogenesis_rate_g_per_h: 2.5, + beta_oxidation_rate_g_per_h: 4.5, + ammonia_clearance_umol_min: 750.0, + albumin_g_dl: 4.1, + clotting_factor_synthesis_pct: 100.0, + bile_acid_synthesis_mg_min: 9.0, + bile_secretion_ml_min: 0.75, + kupffer_activation: 0.25, + acute_phase_response: 0.1, + hepatic_blood_flow_l_min: 1.35, + portal_pressure_mm_hg: 7.0, + hepatic_fat_fraction_pct: 8.0, + insulin_signal: 0.35, + glucagon_signal: 0.45, + cortisol_signal: 0.25, + oxidative_stress_index: 0.18, + state: HepaticState::Postabsorptive, + time_in_state_s: 0.0, + feeding_clock_s: 0.0, + target_meal_interval_s: 4.5 * 3600.0, } } + + fn approach(current: f32, target: f32, rate_per_second: f32, dt_seconds: f32) -> f32 { + let rate = rate_per_second.max(0.0); + if rate == 0.0 || dt_seconds <= 0.0 { + return current; + } + let delta = target - current; + let max_step = rate * dt_seconds; + if delta > max_step { + current + max_step + } else if delta < -max_step { + current - max_step + } else { + target + } + } + + fn simulate_hormone_inputs(&mut self, dt_seconds: f32) { + self.feeding_clock_s += dt_seconds; + if self.feeding_clock_s >= self.target_meal_interval_s { + self.insulin_signal = 0.95; + self.glucagon_signal = 0.2; + self.feeding_clock_s = 0.0; + self.target_meal_interval_s = + (4.0 + 1.5 * (self.hepatic_fat_fraction_pct / 20.0)) * 3600.0; + self.time_in_state_s = 0.0; + self.state = HepaticState::FedAnabolic; + } else { + self.insulin_signal = Self::approach(self.insulin_signal, 0.3, 0.1, dt_seconds); + self.glucagon_signal = Self::approach(self.glucagon_signal, 0.55, 0.08, dt_seconds); + } + self.cortisol_signal = Self::approach( + self.cortisol_signal, + (0.25 + + 0.3 * (self.glucagon_signal - 0.5).max(0.0) + + 0.2 * self.oxidative_stress_index) + .clamp(0.1, 0.9), + 0.05, + dt_seconds, + ); + } + + fn transition_state(&mut self, new_state: HepaticState) { + if self.state != new_state { + self.state = new_state; + self.time_in_state_s = 0.0; + } + } + + fn update_state(&mut self) { + match self.state { + HepaticState::Postabsorptive => { + if self.oxidative_stress_index > 0.6 && self.glycogen_store_g < 40.0 { + self.transition_state(HepaticState::Regenerating); + } else if self.insulin_signal > 0.6 { + self.transition_state(HepaticState::FedAnabolic); + } else if self.glucagon_signal > 0.7 { + self.transition_state(HepaticState::FastingCatabolic); + } else if self.acute_phase_response > 0.4 { + self.transition_state(HepaticState::AcutePhaseResponse); + } + } + HepaticState::FedAnabolic => { + if self.time_in_state_s > 2.0 * 3600.0 && self.glycogen_store_g > 70.0 { + self.transition_state(HepaticState::Postabsorptive); + } + } + HepaticState::FastingCatabolic => { + if self.oxidative_stress_index > 0.7 && self.glycogen_store_g < 30.0 { + self.transition_state(HepaticState::Regenerating); + } else if self.insulin_signal > 0.55 { + self.transition_state(HepaticState::FedAnabolic); + } else if self.time_in_state_s > 12.0 * 3600.0 { + self.transition_state(HepaticState::Postabsorptive); + } + } + HepaticState::AcutePhaseResponse => { + if self.acute_phase_response < 0.2 { + self.transition_state(HepaticState::Postabsorptive); + } else if self.oxidative_stress_index > 0.65 { + self.transition_state(HepaticState::Regenerating); + } + } + HepaticState::Regenerating => { + if self.oxidative_stress_index < 0.3 && self.glycogen_store_g > 60.0 { + self.transition_state(HepaticState::Postabsorptive); + } + } + } + } + fn update_metabolic_fluxes(&mut self, dt_seconds: f32) { + let ( + gluconeogenesis_target, + glycogenolysis_target, + lipogenesis_target, + beta_oxidation_target, + ) = match self.state { + HepaticState::FedAnabolic => (0.8, 2.0, 6.0, 2.0), + HepaticState::Postabsorptive => (1.8, 6.0, 2.5, 4.5), + HepaticState::FastingCatabolic => (2.6, 9.0, 1.2, 6.5), + HepaticState::AcutePhaseResponse => (2.1, 5.0, 1.8, 4.0), + HepaticState::Regenerating => (1.5, 4.0, 3.5, 3.5), + }; + self.gluconeogenesis_rate = Self::approach( + self.gluconeogenesis_rate, + (gluconeogenesis_target + 0.6 * (self.cortisol_signal - 0.3)).clamp(0.4, 3.5), + 0.05, + dt_seconds, + ); + self.glycogenolysis_rate_g_per_h = Self::approach( + self.glycogenolysis_rate_g_per_h, + (glycogenolysis_target + 4.0 * (self.glucagon_signal - self.insulin_signal)) + .clamp(0.0, 14.0), + 0.1, + dt_seconds, + ); + self.lipogenesis_rate_g_per_h = Self::approach( + self.lipogenesis_rate_g_per_h, + (lipogenesis_target + 5.0 * (self.insulin_signal - 0.4).max(0.0)).clamp(0.5, 12.0), + 0.06, + dt_seconds, + ); + self.beta_oxidation_rate_g_per_h = Self::approach( + self.beta_oxidation_rate_g_per_h, + (beta_oxidation_target + 3.5 * (self.glucagon_signal - 0.5).max(0.0)).clamp(1.0, 10.0), + 0.08, + dt_seconds, + ); + let glycogen_change = (self.lipogenesis_rate_g_per_h * 0.2 + - self.glycogenolysis_rate_g_per_h + - self.gluconeogenesis_rate * 0.6) + * dt_seconds + / 3600.0; + self.glycogen_store_g = (self.glycogen_store_g + glycogen_change).clamp(10.0, 140.0); + self.oxidative_stress_index = Self::approach( + self.oxidative_stress_index, + (0.15 + + 0.25 * (self.beta_oxidation_rate_g_per_h - 3.0) / 7.0 + + 0.2 * self.kupffer_activation) + .clamp(0.05, 0.9), + 0.02, + dt_seconds, + ); + self.metabolism = (self.gluconeogenesis_rate / 1.8 + + self.lipogenesis_rate_g_per_h / 3.0 + + self.beta_oxidation_rate_g_per_h / 4.5) + / 3.0; + } + + fn update_bile_and_detox(&mut self, dt_seconds: f32) { + let bile_target = (0.75 + + 0.3 * (self.lipogenesis_rate_g_per_h / 6.0) + + 0.2 * (self.kupffer_activation - 0.3)) + .clamp(0.3, 1.5); + self.bile_secretion_ml_min = + Self::approach(self.bile_secretion_ml_min, bile_target, 0.04, dt_seconds); + self.bile_acid_synthesis_mg_min = Self::approach( + self.bile_acid_synthesis_mg_min, + (9.0 + 4.0 * (self.glucagon_signal - 0.4)).clamp(3.0, 20.0), + 0.05, + dt_seconds, + ); + let detox_target = (100.0 - 15.0 * self.oxidative_stress_index + + 10.0 * (self.insulin_signal - 0.4)) + .clamp(40.0, 110.0); + self.detox = detox_target.round() as u8; + self.enzymes = Self::approach( + self.enzymes, + (1.0 + 0.4 * (self.cortisol_signal - 0.3) + 0.5 * (self.oxidative_stress_index - 0.2)) + .clamp(0.4, 1.8), + 0.03, + dt_seconds, + ); + self.ammonia_clearance_umol_min = Self::approach( + self.ammonia_clearance_umol_min, + (750.0 + 120.0 * (self.metabolism - 1.0) - 200.0 * self.oxidative_stress_index) + .clamp(200.0, 900.0), + 0.2, + dt_seconds, + ); + } + + fn update_hemodynamics(&mut self, dt_seconds: f32) { + let flow_target = (1.35 + + 0.3 * (self.portal_pressure_mm_hg - 7.0) / 4.0 + + 0.25 * (self.metabolism - 1.0)) + .clamp(0.8, 2.0); + self.hepatic_blood_flow_l_min = + Self::approach(self.hepatic_blood_flow_l_min, flow_target, 0.05, dt_seconds); + self.portal_pressure_mm_hg = Self::approach( + self.portal_pressure_mm_hg, + (7.0 + 1.5 * (self.hepatic_fat_fraction_pct / 10.0 - 0.8) + + 0.8 * (self.kupffer_activation - 0.3)) + .clamp(4.0, 16.0), + 0.05, + dt_seconds, + ); + self.kupffer_activation = Self::approach( + self.kupffer_activation, + (0.25 + 0.4 * self.acute_phase_response + 0.2 * self.portal_pressure_mm_hg / 15.0) + .clamp(0.1, 0.95), + 0.04, + dt_seconds, + ); + self.acute_phase_response = Self::approach( + self.acute_phase_response, + (0.1 + 0.6 * (self.oxidative_stress_index - 0.2).max(0.0)).clamp(0.05, 0.9), + 0.03, + dt_seconds, + ); + } + + fn update_proteins(&mut self, dt_seconds: f32) { + self.albumin_g_dl = Self::approach( + self.albumin_g_dl, + (4.2 - 0.4 * self.acute_phase_response + 0.2 * (self.insulin_signal - 0.4)) + .clamp(2.5, 4.8), + 0.01, + dt_seconds, + ); + self.clotting_factor_synthesis_pct = Self::approach( + self.clotting_factor_synthesis_pct, + (100.0 + - 20.0 * self.oxidative_stress_index + - 15.0 * (0.5 - self.albumin_g_dl / 4.0).max(0.0)) + .clamp(40.0, 120.0), + 0.05, + dt_seconds, + ); + } + + fn update_fat_fraction(&mut self, dt_seconds: f32) { + let fat_delta = (self.lipogenesis_rate_g_per_h - self.beta_oxidation_rate_g_per_h) + * dt_seconds + / (24.0 * 3600.0); + self.hepatic_fat_fraction_pct = + (self.hepatic_fat_fraction_pct + fat_delta * 100.0).clamp(2.0, 25.0); + } } impl Organ for Liver { @@ -30,16 +344,29 @@ impl Organ for Liver { fn organ_type(&self) -> OrganType { self.info.kind() } - fn update(&mut self, _dt_seconds: f32) { - self.enzymes = self.enzymes.clamp(0.0, 2.0); + fn update(&mut self, dt_seconds: f32) { + if dt_seconds <= 0.0 { + return; + } + + self.time_in_state_s += dt_seconds; + + self.simulate_hormone_inputs(dt_seconds); + self.update_state(); + self.update_metabolic_fluxes(dt_seconds); + self.update_bile_and_detox(dt_seconds); + self.update_hemodynamics(dt_seconds); + self.update_proteins(dt_seconds); + self.update_fat_fraction(dt_seconds); } fn summary(&self) -> String { format!( - "Liver[id={}, detox={}, k={:.2}, enz={:.2}]", + "Liver[id={}, state={:?}, glycogen={:.0} g, albumin={:.1} g/dL, bile={:.2} ml/min]", self.id(), - self.detox, - self.metabolism, - self.enzymes + self.state, + self.glycogen_store_g, + self.albumin_g_dl, + self.bile_secretion_ml_min ) } fn as_any(&self) -> &dyn core::any::Any { diff --git a/src/organs/lungs.rs b/src/organs/lungs.rs index eb3ffb5..cd91e58 100644 --- a/src/organs/lungs.rs +++ b/src/organs/lungs.rs @@ -1,7 +1,17 @@ use super::{Organ, OrganInfo}; use crate::types::OrganType; -/// Pulmonary model tracking respiratory rate and oxygen saturation. +/// Ventilatory operating mode reflecting dominant chemoreceptor drive. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum VentilatoryState { + Resting, + HypercapnicResponse, + HypoxicResponse, + ExerciseAugmented, + MechanicalDistress, +} + +/// Pulmonary model tracking ventilation mechanics and gas exchange. #[derive(Debug, Clone)] pub struct Lungs { info: OrganInfo, @@ -9,8 +19,45 @@ pub struct Lungs { pub respiratory_rate_bpm: f32, /// Peripheral oxygen saturation percent. pub spo2_pct: f32, - /// Respiratory distress flag reduces SpO2. + /// Flag indicating external distress/vq mismatch triggers. pub distress: bool, + /// Tidal volume (ml). + pub tidal_volume_ml: f32, + /// Minute ventilation (L/min). + pub minute_ventilation_l_min: f32, + /// Dead-space fraction of each breath (0..=0.5). + pub dead_space_fraction: f32, + /// Alveolar oxygen partial pressure (mmHg). + pub alveolar_po2_mm_hg: f32, + /// Alveolar carbon dioxide partial pressure (mmHg). + pub alveolar_pco2_mm_hg: f32, + /// End tidal CO2 (mmHg). + pub end_tidal_co2_mm_hg: f32, + /// Lung compliance (ml/cmH2O). + pub compliance_ml_cm_h2o: f32, + /// Airway resistance (cmH2O·s/L). + pub airway_resistance_cm_h2o_l_s: f32, + /// Respiratory muscle drive (0..=1). + pub muscle_drive: f32, + /// Chemoreceptor drive (0..=1). + pub chemoreceptor_drive: f32, + /// Ventilation/perfusion ratio. + pub ventilation_perfusion_ratio: f32, + /// Shunt fraction (0..=0.4). + pub shunt_fraction: f32, + /// Pulmonary artery pressure (mmHg). + pub pulmonary_artery_pressure_mm_hg: f32, + /// Pulmonary capillary wedge pressure (mmHg). + pub pcwp_mm_hg: f32, + /// Oxygen delivery (ml O2/min). + pub oxygen_delivery_ml_min: f32, + /// CO2 elimination (ml/min). + pub co2_elimination_ml_min: f32, + /// Functional state. + pub state: VentilatoryState, + time_in_state_s: f32, + metabolic_o2_consumption_ml_min: f32, + metabolic_co2_production_ml_min: f32, } impl Lungs { @@ -21,8 +68,233 @@ impl Lungs { respiratory_rate_bpm: 14.0, spo2_pct: 98.0, distress: false, + tidal_volume_ml: 500.0, + minute_ventilation_l_min: 6.5, + dead_space_fraction: 0.28, + alveolar_po2_mm_hg: 100.0, + alveolar_pco2_mm_hg: 38.0, + end_tidal_co2_mm_hg: 36.0, + compliance_ml_cm_h2o: 110.0, + airway_resistance_cm_h2o_l_s: 2.0, + muscle_drive: 0.45, + chemoreceptor_drive: 0.4, + ventilation_perfusion_ratio: 0.96, + shunt_fraction: 0.03, + pulmonary_artery_pressure_mm_hg: 18.0, + pcwp_mm_hg: 9.0, + oxygen_delivery_ml_min: 960.0, + co2_elimination_ml_min: 180.0, + state: VentilatoryState::Resting, + time_in_state_s: 0.0, + metabolic_o2_consumption_ml_min: 250.0, + metabolic_co2_production_ml_min: 200.0, } } + + fn approach(current: f32, target: f32, rate_per_second: f32, dt_seconds: f32) -> f32 { + let rate = rate_per_second.max(0.0); + if rate == 0.0 || dt_seconds <= 0.0 { + return current; + } + let delta = target - current; + let max_step = rate * dt_seconds; + if delta > max_step { + current + max_step + } else if delta < -max_step { + current - max_step + } else { + target + } + } + + fn update_metabolic_demand(&mut self, dt_seconds: f32) { + let exercise_factor = + matches!(self.state, VentilatoryState::ExerciseAugmented) as u8 as f32; + let distress_factor = if self.distress { 0.2 } else { 0.0 }; + let o2_target = 250.0 * (1.0 + 0.8 * exercise_factor + distress_factor); + let co2_target = 200.0 * (1.0 + 0.9 * exercise_factor + distress_factor); + self.metabolic_o2_consumption_ml_min = Self::approach( + self.metabolic_o2_consumption_ml_min, + o2_target, + 0.4, + dt_seconds, + ); + self.metabolic_co2_production_ml_min = Self::approach( + self.metabolic_co2_production_ml_min, + co2_target, + 0.4, + dt_seconds, + ); + } + + fn update_state(&mut self) { + self.state = if self.distress { + VentilatoryState::MechanicalDistress + } else if self.alveolar_pco2_mm_hg > 45.0 { + VentilatoryState::HypercapnicResponse + } else if self.alveolar_po2_mm_hg < 70.0 { + VentilatoryState::HypoxicResponse + } else if self.metabolic_o2_consumption_ml_min > 300.0 { + VentilatoryState::ExerciseAugmented + } else { + VentilatoryState::Resting + }; + } + + fn chemoreceptor_targets(&self) -> (f32, f32) { + match self.state { + VentilatoryState::Resting => (0.45, 0.48), + VentilatoryState::HypercapnicResponse => (0.8, 0.75), + VentilatoryState::HypoxicResponse => (0.85, 0.82), + VentilatoryState::ExerciseAugmented => (0.9, 0.7), + VentilatoryState::MechanicalDistress => (0.95, 0.65), + } + } + + fn update_drives(&mut self, dt_seconds: f32) { + let (chemo_target, muscle_target) = self.chemoreceptor_targets(); + let hypoxia_error = (90.0 - self.alveolar_po2_mm_hg).max(0.0) / 40.0; + let hypercapnia_error = (self.alveolar_pco2_mm_hg - 40.0).max(0.0) / 20.0; + let drive_boost = (hypoxia_error + hypercapnia_error).clamp(0.0, 1.0); + self.chemoreceptor_drive = Self::approach( + self.chemoreceptor_drive, + (chemo_target + 0.6 * drive_boost).clamp(0.2, 1.0), + 0.8, + dt_seconds, + ); + self.muscle_drive = Self::approach( + self.muscle_drive, + (muscle_target + 0.5 * drive_boost).clamp(0.2, 1.0), + 0.6, + dt_seconds, + ); + } + + fn update_mechanics(&mut self, dt_seconds: f32) { + let rate_target = (12.0 + + 18.0 * self.muscle_drive + + 6.0 * (self.chemoreceptor_drive - 0.5).max(0.0) + + 8.0 * matches!(self.state, VentilatoryState::MechanicalDistress) as i32 as f32) + .clamp(8.0, 40.0); + self.respiratory_rate_bpm = + Self::approach(self.respiratory_rate_bpm, rate_target, 1.5, dt_seconds); + let compliance_target = if self.distress { + 65.0 + } else { + 110.0 - 20.0 * (self.muscle_drive - 0.5).max(0.0) + } + .clamp(40.0, 140.0); + self.compliance_ml_cm_h2o = Self::approach( + self.compliance_ml_cm_h2o, + compliance_target, + 0.2, + dt_seconds, + ); + let resistance_target = if self.distress { + 4.5 + } else { + 2.0 + 1.5 * (0.4 - self.compliance_ml_cm_h2o / 150.0).max(0.0) + } + .clamp(1.2, 6.0); + self.airway_resistance_cm_h2o_l_s = Self::approach( + self.airway_resistance_cm_h2o_l_s, + resistance_target, + 0.3, + dt_seconds, + ); + let tidal_target = (450.0 + 160.0 * (self.muscle_drive - 0.4) + - 50.0 * self.airway_resistance_cm_h2o_l_s) + .clamp(250.0, 900.0); + self.tidal_volume_ml = Self::approach(self.tidal_volume_ml, tidal_target, 30.0, dt_seconds); + self.dead_space_fraction = Self::approach( + self.dead_space_fraction, + (0.28 + + 0.15 * (self.airway_resistance_cm_h2o_l_s - 2.0).max(0.0) + + 0.1 * self.shunt_fraction) + .clamp(0.15, 0.5), + 0.2, + dt_seconds, + ); + let alveolar_ventilation = (self.tidal_volume_ml * (1.0 - self.dead_space_fraction) + / 1000.0) + * self.respiratory_rate_bpm; + self.minute_ventilation_l_min = + (self.tidal_volume_ml / 1000.0 * self.respiratory_rate_bpm).clamp(3.0, 25.0); + self.ventilation_perfusion_ratio = Self::approach( + self.ventilation_perfusion_ratio, + (alveolar_ventilation / 5.0).clamp(0.4, 1.4), + 0.3, + dt_seconds, + ); + } + + fn update_gas_exchange(&mut self, dt_seconds: f32) { + let effective_ventilation = + self.minute_ventilation_l_min * (1.0 - self.dead_space_fraction); + let po2_target = (100.0 + 12.0 * (effective_ventilation - 5.5) + - 30.0 * self.shunt_fraction + - 15.0 * (1.0 - self.ventilation_perfusion_ratio)) + .clamp(40.0, 120.0); + let pco2_target = (40.0 - 5.0 * (effective_ventilation - 6.0) + 10.0 * self.shunt_fraction) + .clamp(25.0, 60.0); + self.alveolar_po2_mm_hg = + Self::approach(self.alveolar_po2_mm_hg, po2_target, 0.5, dt_seconds); + self.alveolar_pco2_mm_hg = + Self::approach(self.alveolar_pco2_mm_hg, pco2_target, 0.5, dt_seconds); + self.end_tidal_co2_mm_hg = Self::approach( + self.end_tidal_co2_mm_hg, + self.alveolar_pco2_mm_hg, + 1.2, + dt_seconds, + ); + let spo2_target = (97.0 + 8.0 * (self.alveolar_po2_mm_hg - 90.0) / 40.0 + - 12.0 * self.shunt_fraction + - 5.0 * (self.metabolic_o2_consumption_ml_min - 250.0) / 200.0) + .clamp(70.0, 100.0); + self.spo2_pct = Self::approach(self.spo2_pct, spo2_target, 0.6, dt_seconds); + self.shunt_fraction = Self::approach( + self.shunt_fraction, + (0.03 + + 0.2 * (1.0 - self.ventilation_perfusion_ratio).max(0.0) + + if self.distress { 0.12 } else { 0.0 }) + .clamp(0.0, 0.35), + 0.4, + dt_seconds, + ); + self.oxygen_delivery_ml_min = Self::approach( + self.oxygen_delivery_ml_min, + self.spo2_pct * 10.0, + 2.0, + dt_seconds, + ); + self.co2_elimination_ml_min = Self::approach( + self.co2_elimination_ml_min, + (self.metabolic_co2_production_ml_min + * (self.minute_ventilation_l_min / 6.0).clamp(0.5, 2.0)) + .clamp(80.0, 600.0), + 1.5, + dt_seconds, + ); + } + + fn update_pressures(&mut self, dt_seconds: f32) { + let pap_target = (18.0 + + 8.0 * (self.shunt_fraction - 0.05).max(0.0) + + 4.0 * (self.minute_ventilation_l_min - 6.0) / 6.0) + .clamp(12.0, 35.0); + self.pulmonary_artery_pressure_mm_hg = Self::approach( + self.pulmonary_artery_pressure_mm_hg, + pap_target, + 0.2, + dt_seconds, + ); + self.pcwp_mm_hg = Self::approach( + self.pcwp_mm_hg, + (8.0 + 0.5 * (self.pulmonary_artery_pressure_mm_hg - 18.0)).clamp(5.0, 18.0), + 0.2, + dt_seconds, + ); + } } impl Organ for Lungs { @@ -33,22 +305,26 @@ impl Organ for Lungs { self.info.kind() } fn update(&mut self, dt_seconds: f32) { - let dt = dt_seconds.clamp(0.0, 10.0); - let target_rr = 14.0; - self.respiratory_rate_bpm += 0.1 * (target_rr - self.respiratory_rate_bpm) * (dt / 1.0); - // distress drifts SpO2 downward - if self.distress { - self.spo2_pct -= 0.5 * (dt / 1.0); + if dt_seconds <= 0.0 { + return; } - // keep spo2 in [70, 100] - self.spo2_pct = self.spo2_pct.clamp(70.0, 100.0); + self.time_in_state_s += dt_seconds; + self.update_metabolic_demand(dt_seconds); + self.update_state(); + self.update_drives(dt_seconds); + self.update_mechanics(dt_seconds); + self.update_gas_exchange(dt_seconds); + self.update_pressures(dt_seconds); } fn summary(&self) -> String { format!( - "Lungs[id={}, RR={:.1} bpm, SpO2={:.0}%]", + "Lungs[id={}, state={:?}, RR={:.0}, VT={:.0} ml, SpO2={:.0}%, PaO2~{:.0}]", self.id(), + self.state, self.respiratory_rate_bpm, - self.spo2_pct + self.tidal_volume_ml, + self.spo2_pct, + self.alveolar_po2_mm_hg ) } fn as_any(&self) -> &dyn core::any::Any { diff --git a/src/organs/pancreas.rs b/src/organs/pancreas.rs index ceef908..e06ebf1 100644 --- a/src/organs/pancreas.rs +++ b/src/organs/pancreas.rs @@ -1,20 +1,206 @@ use super::{Organ, OrganInfo}; use crate::types::OrganType; +/// Dominant endocrine/exocrine activity mode of the pancreas. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum PancreaticState { + Basal, + PostprandialAnabolic, + HypoglycemicCounterregulation, + BetaCellExhaustion, +} + #[derive(Debug, Clone)] pub struct Pancreas { info: OrganInfo, - /// Insulin secretion index + /// Insulin secretion index (µU/mL proxy). pub insulin: f32, + /// Glucagon secretion index (pg/mL proxy). + pub glucagon: f32, + /// Somatostatin output (pg/mL proxy). + pub somatostatin: f32, + /// Pancreatic polypeptide level. + pub pancreatic_polypeptide: f32, + /// Enzyme output (kIU/min). + pub digestive_enzyme_output: f32, + /// Bicarbonate secretion (mmol/min). + pub bicarbonate_output_mmol_min: f32, + /// Estimated beta-cell functional mass fraction (0..=1). + pub beta_cell_mass_fraction: f32, + /// Islet stress index (0..=1). + pub islet_stress_index: f32, + /// Acinar secretion flow (ml/min). + pub acinar_flow_ml_min: f32, + /// Ductal pressure (cmH2O). + pub duct_pressure_cm_h2o: f32, + /// Blood glucose sensed by islets (mg/dL). + pub blood_glucose_mg_dl: f32, + /// Incretin stimulus (0..=1). + pub incretin_signal: f32, + /// Autonomic tone (-1 vagal, +1 sympathetic). + pub autonomic_tone: f32, + /// Current pancreas state. + pub state: PancreaticState, + time_in_state_s: f32, + feeding_clock_s: f32, + target_meal_interval_s: f32, } impl Pancreas { pub fn new(id: impl Into) -> Self { Self { info: OrganInfo::new(id, OrganType::Pancreas), - insulin: 1.0, + insulin: 12.0, + glucagon: 60.0, + somatostatin: 20.0, + pancreatic_polypeptide: 120.0, + digestive_enzyme_output: 18.0, + bicarbonate_output_mmol_min: 1.8, + beta_cell_mass_fraction: 0.92, + islet_stress_index: 0.25, + acinar_flow_ml_min: 0.7, + duct_pressure_cm_h2o: 6.0, + blood_glucose_mg_dl: 95.0, + incretin_signal: 0.2, + autonomic_tone: 0.0, + state: PancreaticState::Basal, + time_in_state_s: 0.0, + feeding_clock_s: 0.0, + target_meal_interval_s: 4.2 * 3600.0, } } + + fn approach(current: f32, target: f32, rate_per_second: f32, dt_seconds: f32) -> f32 { + let rate = rate_per_second.max(0.0); + if rate == 0.0 || dt_seconds <= 0.0 { + return current; + } + let delta = target - current; + let max_step = rate * dt_seconds; + if delta > max_step { + current + max_step + } else if delta < -max_step { + current - max_step + } else { + target + } + } + + fn simulate_meals(&mut self, dt_seconds: f32) { + self.feeding_clock_s += dt_seconds; + if self.feeding_clock_s >= self.target_meal_interval_s { + self.blood_glucose_mg_dl = 155.0; + self.incretin_signal = 0.85; + self.autonomic_tone = -0.4; // vagal dominance + self.feeding_clock_s = 0.0; + self.state = PancreaticState::PostprandialAnabolic; + self.time_in_state_s = 0.0; + self.target_meal_interval_s = (3.5 + 1.2 * self.islet_stress_index) * 3600.0; + } else { + self.incretin_signal = Self::approach(self.incretin_signal, 0.15, 0.06, dt_seconds); + self.autonomic_tone = Self::approach(self.autonomic_tone, 0.1, 0.08, dt_seconds); + } + self.blood_glucose_mg_dl = Self::approach( + self.blood_glucose_mg_dl, + 90.0 + 12.0 * (-self.autonomic_tone).max(0.0), + 0.1, + dt_seconds, + ); + } + + fn update_state(&mut self) { + self.state = if self.beta_cell_mass_fraction < 0.6 || self.islet_stress_index > 0.75 { + PancreaticState::BetaCellExhaustion + } else if self.blood_glucose_mg_dl < 70.0 { + PancreaticState::HypoglycemicCounterregulation + } else if self.blood_glucose_mg_dl > 130.0 || self.incretin_signal > 0.5 { + PancreaticState::PostprandialAnabolic + } else { + PancreaticState::Basal + }; + } + + fn update_endocrine(&mut self, dt_seconds: f32) { + let insulin_target = match self.state { + PancreaticState::PostprandialAnabolic => { + 8.0 + 0.6 * (self.blood_glucose_mg_dl - 90.0).max(0.0) + 25.0 * self.incretin_signal + } + PancreaticState::Basal => 10.0 + 0.2 * (self.blood_glucose_mg_dl - 90.0), + PancreaticState::HypoglycemicCounterregulation => 4.0, + PancreaticState::BetaCellExhaustion => 6.0, + }; + self.insulin = Self::approach( + self.insulin, + (insulin_target * self.beta_cell_mass_fraction).clamp(2.0, 80.0), + 0.5, + dt_seconds, + ); + let glucagon_target = match self.state { + PancreaticState::HypoglycemicCounterregulation => 150.0, + PancreaticState::Basal => 70.0, + PancreaticState::PostprandialAnabolic => 40.0, + PancreaticState::BetaCellExhaustion => 110.0, + }; + self.glucagon = Self::approach( + self.glucagon, + (glucagon_target + 20.0 * self.autonomic_tone.max(0.0)).clamp(20.0, 200.0), + 0.4, + dt_seconds, + ); + let somatostatin_target = (20.0 + + 15.0 * (self.incretin_signal - 0.3).max(0.0) + + 0.3 * (self.blood_glucose_mg_dl - 90.0)) + .clamp(10.0, 80.0); + self.somatostatin = Self::approach(self.somatostatin, somatostatin_target, 0.5, dt_seconds); + self.pancreatic_polypeptide = Self::approach( + self.pancreatic_polypeptide, + (100.0 + 80.0 * (-self.autonomic_tone).max(0.0) + 40.0 * self.incretin_signal) + .clamp(60.0, 260.0), + 0.3, + dt_seconds, + ); + self.islet_stress_index = Self::approach( + self.islet_stress_index, + (0.2 + 0.4 * (self.blood_glucose_mg_dl - 100.0).max(0.0) / 80.0 + + 0.3 * (self.autonomic_tone).max(0.0)) + .clamp(0.05, 0.95), + 0.04, + dt_seconds, + ); + self.beta_cell_mass_fraction = (self.beta_cell_mass_fraction + - 0.00002 * dt_seconds * (self.islet_stress_index - 0.3).max(0.0) + + 0.000015 * dt_seconds * (0.5 - self.islet_stress_index).max(0.0)) + .clamp(0.4, 1.05); + } + + fn update_exocrine(&mut self, dt_seconds: f32) { + let enzyme_target = + (15.0 + 25.0 * self.incretin_signal + 10.0 * (-self.autonomic_tone).max(0.0)) + .clamp(5.0, 60.0); + self.digestive_enzyme_output = + Self::approach(self.digestive_enzyme_output, enzyme_target, 0.5, dt_seconds); + let bicarb_target = + (1.5 + 2.5 * self.incretin_signal - 0.5 * self.islet_stress_index).clamp(0.5, 5.0); + self.bicarbonate_output_mmol_min = Self::approach( + self.bicarbonate_output_mmol_min, + bicarb_target, + 0.4, + dt_seconds, + ); + self.acinar_flow_ml_min = Self::approach( + self.acinar_flow_ml_min, + (0.6 + 0.5 * self.incretin_signal + 0.3 * (-self.autonomic_tone).max(0.0)) + .clamp(0.3, 2.0), + 0.4, + dt_seconds, + ); + self.duct_pressure_cm_h2o = Self::approach( + self.duct_pressure_cm_h2o, + (6.0 + 4.0 * (self.acinar_flow_ml_min - 0.7)).clamp(4.0, 18.0), + 0.3, + dt_seconds, + ); + } } impl Organ for Pancreas { @@ -24,9 +210,26 @@ impl Organ for Pancreas { fn organ_type(&self) -> OrganType { self.info.kind() } - fn update(&mut self, _dt_seconds: f32) {} + fn update(&mut self, dt_seconds: f32) { + if dt_seconds <= 0.0 { + return; + } + + self.time_in_state_s += dt_seconds; + self.simulate_meals(dt_seconds); + self.update_state(); + self.update_endocrine(dt_seconds); + self.update_exocrine(dt_seconds); + } fn summary(&self) -> String { - format!("Pancreas[id={}, insulin={:.2}]", self.id(), self.insulin) + format!( + "Pancreas[id={}, state={:?}, insulin={:.1}, glucagon={:.0}, enzymes={:.1} kIU/min]", + self.id(), + self.state, + self.insulin, + self.glucagon, + self.digestive_enzyme_output + ) } fn as_any(&self) -> &dyn core::any::Any { self diff --git a/src/organs/spinal_cord.rs b/src/organs/spinal_cord.rs index e7c3f75..8404dff 100644 --- a/src/organs/spinal_cord.rs +++ b/src/organs/spinal_cord.rs @@ -1,12 +1,45 @@ use super::{Organ, OrganInfo}; use crate::types::OrganType; +/// Functional state of spinal cord circuitry. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum SpinalCordState { + Intact, + Concussed, + Inflammatory, + Ischemic, + NeurogenicShock, +} + #[derive(Debug, Clone)] pub struct SpinalCord { info: OrganInfo, /// 0..=100 nerve signal integrity. pub signal_integrity: u8, pub injury: bool, + /// Ascending conduction velocity (m/s). + pub ascending_conduction_velocity: f32, + /// Descending motor conduction (m/s). + pub descending_conduction_velocity: f32, + /// Segmental reflex gain (dimensionless 0..=2). + pub reflex_gain: f32, + /// Sympathetic preganglionic output (0..=1). + pub sympathetic_outflow: f32, + /// Parasympathetic sacral outflow (0..=1). + pub parasympathetic_outflow: f32, + /// Central pattern generator tone (0..=1). + pub locomotor_cpg_tone: f32, + /// Nociceptive facilitation (0..=1). + pub nociceptive_facilitation: f32, + /// Glial scar formation index (0..=1). + pub glial_scar_index: f32, + /// Cord perfusion pressure (mmHg). + pub cord_perfusion_pressure_mm_hg: f32, + /// Inflammation marker (0..=1). + pub inflammation_index: f32, + /// State of spinal cord. + pub state: SpinalCordState, + time_in_state_s: f32, } impl SpinalCord { @@ -15,8 +48,131 @@ impl SpinalCord { info: OrganInfo::new(id, OrganType::SpinalCord), signal_integrity: 100, injury: false, + ascending_conduction_velocity: 54.0, + descending_conduction_velocity: 60.0, + reflex_gain: 1.0, + sympathetic_outflow: 0.6, + parasympathetic_outflow: 0.55, + locomotor_cpg_tone: 0.65, + nociceptive_facilitation: 0.2, + glial_scar_index: 0.05, + cord_perfusion_pressure_mm_hg: 75.0, + inflammation_index: 0.1, + state: SpinalCordState::Intact, + time_in_state_s: 0.0, } } + + fn approach(current: f32, target: f32, rate_per_second: f32, dt_seconds: f32) -> f32 { + let rate = rate_per_second.max(0.0); + if rate == 0.0 || dt_seconds <= 0.0 { + return current; + } + let delta = target - current; + let max_step = rate * dt_seconds; + if delta > max_step { + current + max_step + } else if delta < -max_step { + current - max_step + } else { + target + } + } + + fn update_state(&mut self) { + self.state = if self.cord_perfusion_pressure_mm_hg < 60.0 { + SpinalCordState::Ischemic + } else if self.injury && self.time_in_state_s < 6.0 * 3600.0 { + SpinalCordState::NeurogenicShock + } else if self.inflammation_index > 0.5 { + SpinalCordState::Inflammatory + } else if self.glial_scar_index > 0.4 { + SpinalCordState::Concussed + } else { + SpinalCordState::Intact + }; + } + + fn update_integrity(&mut self, dt_seconds: f32) { + if self.injury { + let drop = (0.03 * dt_seconds).min(self.signal_integrity as f32); + self.signal_integrity = self.signal_integrity.saturating_sub(drop as u8); + self.glial_scar_index = (self.glial_scar_index + 0.00005 * dt_seconds).clamp(0.0, 1.0); + } else { + self.signal_integrity = (self.signal_integrity + 1).min(100); + self.glial_scar_index = Self::approach(self.glial_scar_index, 0.05, 0.0002, dt_seconds); + } + } + + fn update_conduction(&mut self, dt_seconds: f32) { + let integrity_factor = self.signal_integrity as f32 / 100.0; + let scar_penalty = self.glial_scar_index * 20.0; + self.ascending_conduction_velocity = Self::approach( + self.ascending_conduction_velocity, + (54.0 * integrity_factor - scar_penalty).clamp(10.0, 60.0), + 0.2, + dt_seconds, + ); + self.descending_conduction_velocity = Self::approach( + self.descending_conduction_velocity, + (60.0 * integrity_factor - scar_penalty * 1.2).clamp(15.0, 70.0), + 0.2, + dt_seconds, + ); + } + + fn update_autonomic_outflow(&mut self, dt_seconds: f32) { + let (sym_target, parasym_target, reflex_target, nocice_target) = match self.state { + SpinalCordState::Intact => (0.6, 0.55, 1.0, 0.2), + SpinalCordState::Concussed => (0.5, 0.45, 0.8, 0.35), + SpinalCordState::Inflammatory => (0.65, 0.4, 1.1, 0.6), + SpinalCordState::Ischemic => (0.45, 0.35, 0.6, 0.7), + SpinalCordState::NeurogenicShock => (0.3, 0.25, 0.4, 0.5), + }; + self.sympathetic_outflow = + Self::approach(self.sympathetic_outflow, sym_target, 0.4, dt_seconds); + self.parasympathetic_outflow = Self::approach( + self.parasympathetic_outflow, + parasym_target, + 0.4, + dt_seconds, + ); + self.reflex_gain = Self::approach(self.reflex_gain, reflex_target, 0.3, dt_seconds); + self.nociceptive_facilitation = Self::approach( + self.nociceptive_facilitation, + nocice_target, + 0.3, + dt_seconds, + ); + self.locomotor_cpg_tone = Self::approach( + self.locomotor_cpg_tone, + (0.65 * (self.reflex_gain / 1.0) * (self.descending_conduction_velocity / 60.0)) + .clamp(0.2, 0.9), + 0.2, + dt_seconds, + ); + } + + fn update_perfusion(&mut self, dt_seconds: f32) { + let perfusion_target = (75.0 - 10.0 * (self.sympathetic_outflow - 0.6) + + 6.0 * (self.parasympathetic_outflow - 0.5) + - 15.0 * (1.0 - self.signal_integrity as f32 / 100.0)) + .clamp(40.0, 90.0); + self.cord_perfusion_pressure_mm_hg = Self::approach( + self.cord_perfusion_pressure_mm_hg, + perfusion_target, + 0.2, + dt_seconds, + ); + self.inflammation_index = Self::approach( + self.inflammation_index, + (0.1 + 0.8 * (1.0 - self.cord_perfusion_pressure_mm_hg / 80.0).max(0.0) + + 0.5 * self.glial_scar_index) + .clamp(0.05, 1.0), + 0.02, + dt_seconds, + ); + } } impl Organ for SpinalCord { @@ -26,16 +182,25 @@ impl Organ for SpinalCord { fn organ_type(&self) -> OrganType { self.info.kind() } - fn update(&mut self, _dt_seconds: f32) { - if self.injury { - self.signal_integrity = self.signal_integrity.saturating_sub(1); + fn update(&mut self, dt_seconds: f32) { + if dt_seconds <= 0.0 { + return; } + self.time_in_state_s += dt_seconds; + self.update_integrity(dt_seconds); + self.update_perfusion(dt_seconds); + self.update_state(); + self.update_conduction(dt_seconds); + self.update_autonomic_outflow(dt_seconds); } fn summary(&self) -> String { format!( - "SpinalCord[id={}, integrity={}]", + "SpinalCord[id={}, state={:?}, integrity={}, sym={:.2}, reflex={:.2}]", self.id(), - self.signal_integrity + self.state, + self.signal_integrity, + self.sympathetic_outflow, + self.reflex_gain ) } fn as_any(&self) -> &dyn core::any::Any { diff --git a/src/organs/spleen.rs b/src/organs/spleen.rs index 536a5e4..8f0f756 100644 --- a/src/organs/spleen.rs +++ b/src/organs/spleen.rs @@ -1,11 +1,40 @@ use super::{Organ, OrganInfo}; use crate::types::OrganType; +/// Functional status of the spleen. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum SplenicState { + Homeostatic, + SympatheticContraction, + HyperimmuneActivation, + Sequestration, + Hypofunction, +} + #[derive(Debug, Clone)] pub struct Spleen { info: OrganInfo, - /// Immune activity 0..=100 + /// Immune activity 0..=100. pub immune_activity: u8, + /// Red pulp blood volume (ml). + pub red_pulp_volume_ml: f32, + /// White pulp lymphocyte activation (0..=1). + pub white_pulp_activation: f32, + /// Platelet reservoir (10^9 cells/L contribution). + pub platelet_reservoir: f32, + /// Sympathetic tone (0..=1). + pub sympathetic_tone: f32, + /// Cytokine output (relative units). + pub cytokine_output: f32, + /// Filtered aged erythrocytes (10^6 cells/min). + pub erythrocyte_culling_rate: f32, + /// IgM production (mg/dL). + pub igm_production_mg_dl: f32, + /// Splenic contraction level (0..=1). + pub contraction_fraction: f32, + /// Current spleen state. + pub state: SplenicState, + time_in_state_s: f32, } impl Spleen { @@ -13,8 +42,121 @@ impl Spleen { Self { info: OrganInfo::new(id, OrganType::Spleen), immune_activity: 80, + red_pulp_volume_ml: 180.0, + white_pulp_activation: 0.45, + platelet_reservoir: 70.0, + sympathetic_tone: 0.35, + cytokine_output: 0.2, + erythrocyte_culling_rate: 2.5, + igm_production_mg_dl: 95.0, + contraction_fraction: 0.1, + state: SplenicState::Homeostatic, + time_in_state_s: 0.0, } } + + fn approach(current: f32, target: f32, rate_per_second: f32, dt_seconds: f32) -> f32 { + let rate = rate_per_second.max(0.0); + if rate == 0.0 || dt_seconds <= 0.0 { + return current; + } + let delta = target - current; + let max_step = rate * dt_seconds; + if delta > max_step { + current + max_step + } else if delta < -max_step { + current - max_step + } else { + target + } + } + + fn update_state(&mut self) { + self.state = if self.sympathetic_tone > 0.7 { + SplenicState::SympatheticContraction + } else if self.white_pulp_activation > 0.75 || self.cytokine_output > 0.6 { + SplenicState::HyperimmuneActivation + } else if self.red_pulp_volume_ml > 260.0 { + SplenicState::Sequestration + } else if self.white_pulp_activation < 0.2 { + SplenicState::Hypofunction + } else { + SplenicState::Homeostatic + }; + } + + fn update_sympathetic_tone(&mut self, dt_seconds: f32) { + self.sympathetic_tone = Self::approach( + self.sympathetic_tone, + (0.3 + 0.5 * (self.contraction_fraction - 0.3).max(0.0)).clamp(0.2, 0.95), + 0.3, + dt_seconds, + ); + } + + fn update_contraction(&mut self, dt_seconds: f32) { + let contraction_target = match self.state { + SplenicState::SympatheticContraction => 0.85, + SplenicState::HyperimmuneActivation => 0.45, + SplenicState::Sequestration => 0.15, + SplenicState::Hypofunction => 0.1, + SplenicState::Homeostatic => 0.3, + }; + self.contraction_fraction = Self::approach( + self.contraction_fraction, + contraction_target, + 0.4, + dt_seconds, + ); + let volume_target = (180.0 + 120.0 * (0.4 - self.contraction_fraction)).clamp(80.0, 320.0); + self.red_pulp_volume_ml = + Self::approach(self.red_pulp_volume_ml, volume_target, 0.8, dt_seconds); + self.platelet_reservoir = Self::approach( + self.platelet_reservoir, + (70.0 + 40.0 * (self.red_pulp_volume_ml - 180.0) / 120.0).clamp(20.0, 160.0), + 0.6, + dt_seconds, + ); + } + + fn update_immune_activity(&mut self, dt_seconds: f32) { + let activation_target = match self.state { + SplenicState::HyperimmuneActivation => 0.85, + SplenicState::Sequestration => 0.55, + SplenicState::Hypofunction => 0.18, + SplenicState::SympatheticContraction => 0.4, + SplenicState::Homeostatic => 0.45, + }; + self.white_pulp_activation = Self::approach( + self.white_pulp_activation, + activation_target, + 0.3, + dt_seconds, + ); + self.immune_activity = ((self.white_pulp_activation * 120.0) + + (self.cytokine_output * 40.0)) + .clamp(10.0, 160.0) as u8; + self.cytokine_output = Self::approach( + self.cytokine_output, + (0.2 + 0.8 * (self.white_pulp_activation - 0.3).max(0.0)).clamp(0.05, 1.2), + 0.1, + dt_seconds, + ); + self.erythrocyte_culling_rate = Self::approach( + self.erythrocyte_culling_rate, + (2.0 + 1.5 * (self.red_pulp_volume_ml - 180.0) / 100.0 + 1.2 * self.cytokine_output) + .clamp(0.5, 8.0), + 0.2, + dt_seconds, + ); + self.igm_production_mg_dl = Self::approach( + self.igm_production_mg_dl, + (90.0 + 60.0 * self.white_pulp_activation - 20.0 * self.sympathetic_tone) + .clamp(30.0, 220.0), + 0.4, + dt_seconds, + ); + } } impl Organ for Spleen { @@ -24,9 +166,25 @@ impl Organ for Spleen { fn organ_type(&self) -> OrganType { self.info.kind() } - fn update(&mut self, _dt_seconds: f32) {} + fn update(&mut self, dt_seconds: f32) { + if dt_seconds <= 0.0 { + return; + } + self.time_in_state_s += dt_seconds; + self.update_state(); + self.update_contraction(dt_seconds); + self.update_sympathetic_tone(dt_seconds); + self.update_immune_activity(dt_seconds); + } fn summary(&self) -> String { - format!("Spleen[id={}, immune={}]", self.id(), self.immune_activity) + format!( + "Spleen[id={}, state={:?}, immune={}, redpulp={:.0} ml, platelets={:.0}]", + self.id(), + self.state, + self.immune_activity, + self.red_pulp_volume_ml, + self.platelet_reservoir + ) } fn as_any(&self) -> &dyn core::any::Any { self diff --git a/src/organs/stomach.rs b/src/organs/stomach.rs index 2434213..3299dbe 100644 --- a/src/organs/stomach.rs +++ b/src/organs/stomach.rs @@ -1,11 +1,52 @@ use super::{Organ, OrganInfo}; use crate::types::OrganType; +/// Gastric functional phase. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum GastricPhase { + Fasting, + Cephalic, + Gastric, + Intestinal, + DelayedEmptying, +} + #[derive(Debug, Clone)] pub struct Stomach { info: OrganInfo, - /// Acid level 0..=100 + /// Acid level 0..=100 (higher = more secretion). pub acid_level: u8, + /// Gastric lumen pH. + pub ph: f32, + /// Current gastric volume (ml). + pub volume_ml: f32, + /// Gastric motility index (0..=1). + pub motility_index: f32, + /// Antral pump strength (0..=1). + pub antral_pump_strength: f32, + /// Gastric emptying rate (ml/min). + pub emptying_rate_ml_min: f32, + /// Ghrelin level (pg/mL proxy). + pub ghrelin: f32, + /// Gastrin level (pg/mL proxy). + pub gastrin: f32, + /// Histamine release (relative units). + pub histamine: f32, + /// Somatostatin brake (relative units). + pub somatostatin: f32, + /// Protective mucus production (g/hour). + pub mucus_production_g_per_h: f32, + /// Intrinsic factor secretion (relative units). + pub intrinsic_factor: f32, + /// Vagal tone (0..=1). + pub vagal_tone: f32, + /// Gastric phase. + pub phase: GastricPhase, + /// Pending meal caloric load (kcal). + pub nutrient_load_kcal: f32, + time_in_phase_s: f32, + fasting_clock_s: f32, + target_meal_interval_s: f32, } impl Stomach { @@ -13,8 +54,197 @@ impl Stomach { Self { info: OrganInfo::new(id, OrganType::Stomach), acid_level: 50, + ph: 2.2, + volume_ml: 120.0, + motility_index: 0.35, + antral_pump_strength: 0.3, + emptying_rate_ml_min: 1.5, + ghrelin: 950.0, + gastrin: 80.0, + histamine: 0.4, + somatostatin: 0.3, + mucus_production_g_per_h: 15.0, + intrinsic_factor: 0.6, + vagal_tone: 0.4, + phase: GastricPhase::Fasting, + nutrient_load_kcal: 60.0, + time_in_phase_s: 0.0, + fasting_clock_s: 0.0, + target_meal_interval_s: 4.5 * 3600.0, } } + + fn approach(current: f32, target: f32, rate_per_second: f32, dt_seconds: f32) -> f32 { + let rate = rate_per_second.max(0.0); + if rate == 0.0 || dt_seconds <= 0.0 { + return current; + } + let delta = target - current; + let max_step = rate * dt_seconds; + if delta > max_step { + current + max_step + } else if delta < -max_step { + current - max_step + } else { + target + } + } + + fn simulate_meals(&mut self, dt_seconds: f32) { + self.fasting_clock_s += dt_seconds; + if self.fasting_clock_s >= self.target_meal_interval_s { + self.phase = GastricPhase::Cephalic; + self.time_in_phase_s = 0.0; + self.vagal_tone = 0.85; + self.ghrelin = 600.0; + self.gastrin = 160.0; + self.nutrient_load_kcal = 650.0; + self.volume_ml = (self.volume_ml + 450.0).clamp(80.0, 1600.0); + self.target_meal_interval_s = + (4.0 + 1.0 * (self.mucus_production_g_per_h / 15.0)) * 3600.0; + self.fasting_clock_s = 0.0; + } else { + self.vagal_tone = Self::approach(self.vagal_tone, 0.35, 0.04, dt_seconds); + self.ghrelin = Self::approach(self.ghrelin, 1200.0, 1.0, dt_seconds); + } + } + + fn update_phase(&mut self) { + self.phase = match self.phase { + GastricPhase::Cephalic => { + if self.time_in_phase_s > 300.0 { + GastricPhase::Gastric + } else { + GastricPhase::Cephalic + } + } + GastricPhase::Gastric => { + if self.volume_ml < 200.0 { + GastricPhase::Intestinal + } else { + GastricPhase::Gastric + } + } + GastricPhase::Intestinal => { + if self.nutrient_load_kcal < 80.0 { + GastricPhase::Fasting + } else if self.emptying_rate_ml_min < 1.0 { + GastricPhase::DelayedEmptying + } else { + GastricPhase::Intestinal + } + } + GastricPhase::DelayedEmptying => { + if self.emptying_rate_ml_min > 1.5 { + GastricPhase::Fasting + } else { + GastricPhase::DelayedEmptying + } + } + GastricPhase::Fasting => { + if self.nutrient_load_kcal > 120.0 { + GastricPhase::Cephalic + } else { + GastricPhase::Fasting + } + } + }; + } + + fn update_secretions(&mut self, dt_seconds: f32) { + let gastrin_target = match self.phase { + GastricPhase::Cephalic => 180.0, + GastricPhase::Gastric => 220.0, + GastricPhase::Intestinal => 120.0, + GastricPhase::DelayedEmptying => 160.0, + GastricPhase::Fasting => 60.0, + }; + self.gastrin = Self::approach( + self.gastrin, + (gastrin_target + 0.5 * (self.volume_ml - 250.0).max(0.0)).clamp(40.0, 320.0), + 0.5, + dt_seconds, + ); + self.histamine = Self::approach( + self.histamine, + (0.3 + 0.004 * self.gastrin + 0.2 * (self.vagal_tone - 0.4).max(0.0)).clamp(0.1, 2.0), + 0.3, + dt_seconds, + ); + self.somatostatin = Self::approach( + self.somatostatin, + (0.25 + + 0.2 * (self.ph - 2.0).max(0.0) + + 0.3 * (self.phase == GastricPhase::Intestinal) as i32 as f32) + .clamp(0.1, 2.0), + 0.4, + dt_seconds, + ); + let acid_drive = + (self.gastrin / 200.0 + self.histamine - self.somatostatin).clamp(0.0, 2.0); + let acid_numeric = (50.0 + 35.0 * acid_drive).clamp(10.0, 100.0); + self.acid_level = acid_numeric.round() as u8; + self.ph = Self::approach( + self.ph, + (7.0 - 0.045 * self.acid_level as f32 + 0.4 * (self.volume_ml / 500.0)).clamp(1.2, 6.5), + 0.6, + dt_seconds, + ); + self.mucus_production_g_per_h = Self::approach( + self.mucus_production_g_per_h, + (15.0 + + 6.0 * (self.acid_level as f32 / 60.0) + + 4.0 * (self.somatostatin - 0.3).max(0.0)) + .clamp(8.0, 40.0), + 0.2, + dt_seconds, + ); + self.intrinsic_factor = Self::approach( + self.intrinsic_factor, + (0.6 + 0.4 * (self.acid_level as f32 / 80.0)).clamp(0.2, 1.2), + 0.2, + dt_seconds, + ); + } + + fn update_motility(&mut self, dt_seconds: f32) { + let motility_target = match self.phase { + GastricPhase::Cephalic => 0.4, + GastricPhase::Gastric => 0.75, + GastricPhase::Intestinal => 0.6, + GastricPhase::DelayedEmptying => 0.35, + GastricPhase::Fasting => 0.3, + }; + self.motility_index = Self::approach(self.motility_index, motility_target, 0.5, dt_seconds); + self.antral_pump_strength = Self::approach( + self.antral_pump_strength, + (0.3 + 0.5 * self.motility_index + 0.3 * self.vagal_tone).clamp(0.2, 0.95), + 0.5, + dt_seconds, + ); + self.emptying_rate_ml_min = Self::approach( + self.emptying_rate_ml_min, + (1.5 + 3.5 * self.antral_pump_strength + - 1.0 * (self.ph - 3.0).max(0.0) + - 0.5 * (self.nutrient_load_kcal / 300.0)) + .clamp(0.2, 9.0), + 0.4, + dt_seconds, + ); + } + + fn update_volume(&mut self, dt_seconds: f32) { + let emptied = self.emptying_rate_ml_min * dt_seconds / 60.0; + let metabolic_use = (self.nutrient_load_kcal * 0.3) * dt_seconds / 3600.0; + self.volume_ml = (self.volume_ml - emptied).clamp(30.0, 1800.0); + self.nutrient_load_kcal = (self.nutrient_load_kcal - metabolic_use).max(0.0); + self.ghrelin = Self::approach( + self.ghrelin, + (1200.0 - 0.8 * self.volume_ml).clamp(200.0, 1400.0), + 0.4, + dt_seconds, + ); + } } impl Organ for Stomach { @@ -24,9 +254,26 @@ impl Organ for Stomach { fn organ_type(&self) -> OrganType { self.info.kind() } - fn update(&mut self, _dt_seconds: f32) {} + fn update(&mut self, dt_seconds: f32) { + if dt_seconds <= 0.0 { + return; + } + self.time_in_phase_s += dt_seconds; + self.simulate_meals(dt_seconds); + self.update_phase(); + self.update_secretions(dt_seconds); + self.update_motility(dt_seconds); + self.update_volume(dt_seconds); + } fn summary(&self) -> String { - format!("Stomach[id={}, acid={}]", self.id(), self.acid_level) + format!( + "Stomach[id={}, phase={:?}, vol={:.0} ml, pH={:.1}, acid={}]", + self.id(), + self.phase, + self.volume_ml, + self.ph, + self.acid_level + ) } fn as_any(&self) -> &dyn core::any::Any { self diff --git a/src/patient.rs b/src/patient.rs index 3d12e99..316a0fb 100644 --- a/src/patient.rs +++ b/src/patient.rs @@ -1,8 +1,11 @@ //! Patient type holding organs and core physiology snapshots. use crate::error::MedicalError; -use crate::organs::{Heart, Lungs, Organ}; -use crate::types::{Blood, BloodPressure, OrganType}; +use crate::organs::{ + Bladder, Brain, Gallbladder, Heart, IntestinalPhase, Intestines, Kidneys, Liver, Lungs, Organ, + Pancreas, SpinalCord, Spleen, SplenicState, Stomach, +}; +use crate::organs::intestines::IntestinalPhase;\r\nuse crate::organs::spleen::SplenicState;\r\n\r\nuse crate::types::{Blood, BloodPressure, OrganType}; /// Patient container and simulation entry. #[derive(Debug)] @@ -15,6 +18,68 @@ pub struct Patient { pub blood_pressure: BloodPressure, } +#[derive(Clone, Copy)] +struct HeartSignals { + systolic: f32, + diastolic: f32, + map: f32, + cardiac_output: f32, + heart_rate: f32, +} + +#[derive(Clone, Copy)] +struct LungSignals { + spo2_pct: f32, + alveolar_po2_mm_hg: f32, + alveolar_pco2_mm_hg: f32, + minute_ventilation_l_min: f32, + oxygen_delivery_ml_min: f32, + shunt_fraction: f32, +} + +#[derive(Clone, Copy)] +struct LiverSignals { + bile_secretion_ml_min: f32, +} + +#[derive(Clone, Copy)] +struct PancreasSignals { + insulin: f32, + glucagon: f32, + incretin_signal: f32, + somatostatin: f32, +} + +#[derive(Clone, Copy)] +struct IntestineSignals { + phase: IntestinalPhase, + nutrient_energy_kcal: f32, +} + +#[derive(Clone, Copy)] +struct GallbladderSignals { + bile_acid_concentration_mmol_l: f32, +} + +#[derive(Clone, Copy)] +struct SpinalSignals { + sympathetic_outflow: f32, + parasympathetic_outflow: f32, +} + +#[derive(Clone, Copy)] +struct KidneySignals { + urine_flow_ml_min: f32, +} + +#[derive(Clone, Copy)] +struct SpleenSignals { + immune_activity: u8, + red_pulp_volume_ml: f32, + platelet_reservoir: f32, + state: SplenicState, +} + impl Patient { /// Construct a new patient with validated id. pub fn new(id: impl Into) -> crate::Result { @@ -71,7 +136,7 @@ impl Patient { self.with_heart(12) } - /// Initialize a patient with a heart with `leads`. + /// Initialize a patient with a heart with leads. pub fn with_heart(mut self, leads: u8) -> Self { let id = format!("{}-heart", self.id); self.add_organ(Heart::new(id, leads)); @@ -98,31 +163,31 @@ impl Patient { } OrganType::Brain => { let id = format!("{}-brain", self.id); - self.add_organ(crate::organs::Brain::new(id)); + self.add_organ(Brain::new(id)); } OrganType::SpinalCord => { let id = format!("{}-sc", self.id); - self.add_organ(crate::organs::SpinalCord::new(id)); + self.add_organ(SpinalCord::new(id)); } OrganType::Stomach => { let id = format!("{}-stomach", self.id); - self.add_organ(crate::organs::Stomach::new(id)); + self.add_organ(Stomach::new(id)); } OrganType::Liver => { let id = format!("{}-liver", self.id); - self.add_organ(crate::organs::Liver::new(id)); + self.add_organ(Liver::new(id)); } OrganType::Gallbladder => { let id = format!("{}-gb", self.id); - self.add_organ(crate::organs::Gallbladder::new(id)); + self.add_organ(Gallbladder::new(id)); } OrganType::Pancreas => { let id = format!("{}-pancreas", self.id); - self.add_organ(crate::organs::Pancreas::new(id)); + self.add_organ(Pancreas::new(id)); } OrganType::Intestines => { let id = format!("{}-intestines", self.id); - self.add_organ(crate::organs::Intestines::new(id)); + self.add_organ(Intestines::new(id)); } OrganType::Esophagus => { let id = format!("{}-eso", self.id); @@ -130,42 +195,351 @@ impl Patient { } OrganType::Kidneys => { let id = format!("{}-kidneys", self.id); - self.add_organ(crate::organs::Kidneys::new(id)); + self.add_organ(Kidneys::new(id)); } OrganType::Bladder => { let id = format!("{}-bladder", self.id); - self.add_organ(crate::organs::Bladder::new(id)); + self.add_organ(Bladder::new(id)); } OrganType::Spleen => { let id = format!("{}-spleen", self.id); - self.add_organ(crate::organs::Spleen::new(id)); + self.add_organ(Spleen::new(id)); } } self } - /// Advance simulation by `dt_seconds`. + /// Advance simulation by dt_seconds. pub fn update(&mut self, dt_seconds: f32) { + if dt_seconds <= 0.0 { + return; + } + + let heart_signals = self.find_organ_typed::().map(|h| { + let systolic = h.arterial_bp.systolic as f32; + let diastolic = h.arterial_bp.diastolic as f32; + HeartSignals { + systolic, + diastolic, + map: diastolic + (systolic - diastolic) / 3.0, + cardiac_output: h.cardiac_output_l_min, + heart_rate: h.heart_rate_bpm, + } + }); + + let lungs_signals = self.find_organ_typed::().map(|l| LungSignals { + spo2_pct: l.spo2_pct, + alveolar_po2_mm_hg: l.alveolar_po2_mm_hg, + alveolar_pco2_mm_hg: l.alveolar_pco2_mm_hg, + minute_ventilation_l_min: l.minute_ventilation_l_min, + oxygen_delivery_ml_min: l.oxygen_delivery_ml_min, + shunt_fraction: l.shunt_fraction, + }); + + let liver_signals = self.find_organ_typed::().map(|l| LiverSignals { + bile_secretion_ml_min: l.bile_secretion_ml_min, + }); + + let intestine_signals = self + .find_organ_typed::() + .map(|i| IntestineSignals { + phase: i.phase, + nutrient_energy_kcal: i.nutrient_energy_kcal, + }); + + let gallbladder_signals = + self.find_organ_typed::() + .map(|g| GallbladderSignals { + bile_acid_concentration_mmol_l: g.bile_acid_concentration_mmol_l, + }); + + let kidney_signals = self.find_organ_typed::().map(|k| KidneySignals { + urine_flow_ml_min: k.urine_flow_ml_min, + }); + + let spinal_signals = self + .find_organ_typed::() + .map(|s| SpinalSignals { + sympathetic_outflow: s.sympathetic_outflow, + parasympathetic_outflow: s.parasympathetic_outflow, + }); + + if let Some(pancreas) = self.find_organ_typed_mut::() { + pancreas.blood_glucose_mg_dl = self.blood.glucose_mg_dl; + if let Some(intestines) = intestine_signals { + let incretin_target = (intestines.nutrient_energy_kcal / 400.0).clamp(0.05, 1.0); + pancreas.incretin_signal = + Self::relax_value(pancreas.incretin_signal, incretin_target, dt_seconds, 90.0); + } + if let Some(spinal) = spinal_signals { + let tone_target = + (spinal.sympathetic_outflow - spinal.parasympathetic_outflow).clamp(-1.0, 1.0); + pancreas.autonomic_tone = + Self::relax_value(pancreas.autonomic_tone, tone_target, dt_seconds, 120.0); + } + } + + if let Some(brain) = self.find_organ_typed_mut::() { + if let Some(lungs) = lungs_signals { + brain.oxygenation_saturation = (lungs.spo2_pct / 100.0).clamp(0.8, 1.0); + let drive_target = + (0.55 + (lungs.alveolar_pco2_mm_hg - 38.0) / 40.0).clamp(0.2, 1.0); + brain.brainstem_autonomic_drive = Self::relax_value( + brain.brainstem_autonomic_drive, + drive_target, + dt_seconds, + 20.0, + ); + } + if let Some(heart) = heart_signals { + let cpp_target = (heart.map - brain.intracranial_pressure_mm_hg).clamp(40.0, 110.0); + brain.cerebral_perfusion_pressure_mm_hg = Self::relax_value( + brain.cerebral_perfusion_pressure_mm_hg, + cpp_target, + dt_seconds, + 15.0, + ); + } + } + + if let Some(kidneys) = self.find_organ_typed_mut::() { + let osm_target = 285.0 + (self.blood.glucose_mg_dl - 95.0) * 0.06; + kidneys.serum_osmolality_mosm = + Self::relax_value(kidneys.serum_osmolality_mosm, osm_target, dt_seconds, 120.0); + if let Some(heart) = heart_signals { + let plasma_target = (3.0 + 0.22 * (heart.cardiac_output - 5.0)).clamp(2.4, 3.6); + kidneys.plasma_volume_l = + Self::relax_value(kidneys.plasma_volume_l, plasma_target, dt_seconds, 180.0); + } + } + + if let Some(gallbladder) = self.find_organ_typed_mut::() { + if let Some(liver) = liver_signals { + let inflow_target = (liver.bile_secretion_ml_min * 0.8).clamp(0.05, 2.4); + gallbladder.hepatic_bile_flow_ml_per_min = Self::relax_value( + gallbladder.hepatic_bile_flow_ml_per_min, + inflow_target, + dt_seconds, + 80.0, + ); + } + if let Some(intestines) = intestine_signals { + if matches!( + intestines.phase, + IntestinalPhase::IlealBrake | IntestinalPhase::Dysmotility + ) { + gallbladder.sphincter_of_oddi_tone = (gallbladder.sphincter_of_oddi_tone + + 0.05 * dt_seconds / 60.0) + .clamp(0.2, 0.95); + } + } + } + + if let Some(intestines) = self.find_organ_typed_mut::() { + if let Some(gall) = gallbladder_signals { + let recirc_target = + (0.82 + (gall.bile_acid_concentration_mmol_l - 60.0) / 320.0).clamp(0.5, 0.98); + intestines.bile_acid_recirculation_fraction = Self::relax_value( + intestines.bile_acid_recirculation_fraction, + recirc_target, + dt_seconds, + 240.0, + ); + } + } + + if let Some(spinal) = self.find_organ_typed_mut::() { + if let Some(heart) = heart_signals { + let perfusion_target = (heart.map - 8.0).clamp(45.0, 90.0); + spinal.cord_perfusion_pressure_mm_hg = Self::relax_value( + spinal.cord_perfusion_pressure_mm_hg, + perfusion_target, + dt_seconds, + 120.0, + ); + } + } + + if let Some(bladder) = self.find_organ_typed_mut::() { + if let Some(kidney) = kidney_signals { + let fill_target = (kidney.urine_flow_ml_min * 60.0).clamp(5.0, 180.0); + bladder.filling_rate_ml_per_min = Self::relax_value( + bladder.filling_rate_ml_per_min, + fill_target, + dt_seconds, + 90.0, + ); + } + } + for organ in &mut self.organs { organ.update(dt_seconds); } - // Simple inter-organ signaling: low SpO2 nudges heart rate higher. - if let Some(spo2) = self.find_organ_typed::().map(|l| l.spo2_pct) { - if let Some(heart) = self.find_organ_typed_mut::() { - let target = if spo2 < 92.0 { 90.0 } else { 70.0 }; - let diff = target - heart.heart_rate_bpm; - heart.heart_rate_bpm += 0.05 * diff; + + if let Some(heart) = self.find_organ_typed::() { + self.blood_pressure = heart.arterial_bp; + } + if let Some(lungs) = self.find_organ_typed::() { + self.blood.spo2_pct = lungs.spo2_pct; + } + + let mut glucose = self.blood.glucose_mg_dl; + if let Some(liver) = self.find_organ_typed::() { + let hepatic_balance = liver.gluconeogenesis_rate * 24.0 + + liver.glycogenolysis_rate_g_per_h * 6.0 + - liver.lipogenesis_rate_g_per_h * 4.0 + - liver.beta_oxidation_rate_g_per_h * 2.5 + - (liver.insulin_signal - liver.glucagon_signal) * 30.0; + glucose += hepatic_balance * (dt_seconds / 3600.0); + } + if let Some(pancreas) = self.find_organ_typed::() { + let hormonal_balance = pancreas.glucagon * 0.05 - pancreas.insulin * 0.08; + glucose += hormonal_balance * (dt_seconds / 60.0); + } + glucose = glucose.clamp(60.0, 220.0); + self.blood.glucose_mg_dl = glucose; + + let kidneys_after = self + .find_organ_typed::() + .map(|k| (k.erythropoietin_iu_per_day, k.urine_flow_ml_min)); + + if let Some((epo, _)) = kidneys_after { + let hgb_target = 14.0 + (epo - 18.0) / 80.0; + self.blood.hemoglobin_g_dl = Self::relax_value( + self.blood.hemoglobin_g_dl, + hgb_target.clamp(9.0, 18.0), + dt_seconds, + 600.0, + ); + } + + let spleen_after = self.find_organ_typed::().map(|s| SpleenSignals { + immune_activity: s.immune_activity, + red_pulp_volume_ml: s.red_pulp_volume_ml, + platelet_reservoir: s.platelet_reservoir, + state: s.state, + }); + + if let Some(spleen) = spleen_after { + let state_penalty = match spleen.state { + SplenicState::SympatheticContraction => -2.0, + SplenicState::HyperimmuneActivation => 2.5, + SplenicState::Sequestration => 3.0, + SplenicState::Hypofunction => -1.5, + SplenicState::Homeostatic => 0.0, + }; + let hematocrit_target = 42.0 + - (spleen.red_pulp_volume_ml - 180.0) / 8.0 + - (spleen.platelet_reservoir - 70.0) / 30.0 + + state_penalty; + self.blood.hematocrit_pct = Self::relax_value( + self.blood.hematocrit_pct, + hematocrit_target.clamp(30.0, 55.0), + dt_seconds, + 600.0, + ); + } + + if let Some(pancreas) = self.find_organ_typed_mut::() { + pancreas.blood_glucose_mg_dl = self.blood.glucose_mg_dl; + } + + let pancreas_after = self + .find_organ_typed::() + .map(|p| PancreasSignals { + insulin: p.insulin, + glucagon: p.glucagon, + incretin_signal: p.incretin_signal, + somatostatin: p.somatostatin, + }); + + if let Some(pancreas) = pancreas_after { + if let Some(liver) = self.find_organ_typed_mut::() { + let insulin_target = (pancreas.insulin / 60.0).clamp(0.1, 1.0); + let glucagon_target = (pancreas.glucagon / 120.0).clamp(0.1, 1.2); + liver.insulin_signal = + Self::relax_value(liver.insulin_signal, insulin_target, dt_seconds, 240.0); + liver.glucagon_signal = + Self::relax_value(liver.glucagon_signal, glucagon_target, dt_seconds, 240.0); } } - // Kidneys produce urine into bladder - let produced_opt = self - .find_organ_typed::() - .map(|kidneys| (kidneys.gfr * (dt_seconds / 60.0)).max(0.0) * 0.5); // ml - if let (Some(produced), Some(bladder)) = ( - produced_opt, - self.find_organ_typed_mut::(), - ) { - bladder.volume_ml += produced; + + if let Some(liver) = self.find_organ_typed::() { + if let Some(gallbladder) = self.find_organ_typed_mut::() { + let inflow_target = (liver.bile_secretion_ml_min * 0.8).clamp(0.05, 2.4); + gallbladder.hepatic_bile_flow_ml_per_min = Self::relax_value( + gallbladder.hepatic_bile_flow_ml_per_min, + inflow_target, + dt_seconds, + 80.0, + ); + } + } + + if let Some(stomach) = self.find_organ_typed::() { + let delivered_ml = stomach.emptying_rate_ml_min * dt_seconds / 60.0; + let delivered_kcal = (delivered_ml * 0.8).min(stomach.nutrient_load_kcal); + if delivered_kcal > 0.0 { + if let Some(intestines) = self.find_organ_typed_mut::() { + intestines.nutrient_energy_kcal += delivered_kcal; + } + } + } + + if let Some((_, urine_flow)) = kidneys_after { + let produced = (urine_flow * dt_seconds / 60.0).max(0.0); + if produced > 0.0 { + if let Some(bladder) = self.find_organ_typed_mut::() { + bladder.volume_ml += produced; + } + } + } + + let brain_after = self + .find_organ_typed::() + .map(|b| (b.brainstem_autonomic_drive, b.autonomic_variability)); + let spinal_after = self + .find_organ_typed::() + .map(|s| (s.sympathetic_outflow, s.parasympathetic_outflow)); + if let Some(heart) = self.find_organ_typed_mut::() { + if let Some((brain_drive, brain_sympathetic)) = brain_after { + let mut tone_target = (brain_drive - 0.5) * 1.2 + (brain_sympathetic - 0.5) * 0.8; + if let Some((sym, para)) = spinal_after { + tone_target += (sym - para) * 0.6; + } + heart.autonomic_tone = Self::relax_value( + heart.autonomic_tone, + tone_target.clamp(-1.0, 1.0), + dt_seconds, + 40.0, + ) + .clamp(-1.0, 1.0); + } + } + + if let Some(pancreas) = pancreas_after { + if let Some(stomach) = self.find_organ_typed_mut::() { + let acid_current = stomach.acid_level as f32; + let acid_target = (acid_current - pancreas.somatostatin * 5.0).clamp(10.0, 100.0); + let acid_next = Self::relax_value(acid_current, acid_target, dt_seconds, 120.0); + stomach.acid_level = acid_next.round().clamp(0.0, 100.0) as u8; + stomach.vagal_tone = Self::relax_value( + stomach.vagal_tone, + (0.4 + pancreas.incretin_signal * 0.2).clamp(0.2, 0.9), + dt_seconds, + 120.0, + ); + } + } + } + + #[inline] + fn relax_value(current: f32, target: f32, dt_seconds: f32, time_constant: f32) -> f32 { + if time_constant <= 0.0 { + target + } else { + let alpha = (dt_seconds / time_constant).clamp(0.0, 1.0); + current + (target - current) * alpha } } @@ -210,7 +584,7 @@ fn is_valid_id(id: &str) -> bool { #[cfg(test)] mod tests { use super::*; - use crate::types::OrganType; + use crate::organs::intestines::IntestinalPhase;\r\nuse crate::organs::spleen::SplenicState;\r\n\r\nuse crate::types::OrganType; #[test] fn patient_lifecycle() { @@ -229,3 +603,5 @@ mod tests { assert!(Patient::new("bad id").is_err()); } } + +