organpool 0.2.0

Organ physics substrate — cardiac and respiratory simulation with autonomic modulation, respiratory sinus arrhythmia, integer-only ion channel dynamics, and real-time vital sign diagnostics
Documentation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
//! Respiratory cycle physics — central pattern generator, gas exchange, configuration.
//!
//! The respiratory system models the pre-Botzinger complex as an autonomous
//! oscillator, analogous to the SA node in the cardiac system. A drive ramp
//! accumulates during the end-expiratory pause; when it reaches threshold,
//! inspiration begins. CO2 chemoreceptor input lowers the threshold (more
//! CO2 = breathe sooner). NE increases the ramp rate, ACh decreases it.
//!
//! Physics use wall-clock time — `Instant` timestamps and `Duration` for
//! all timing. Integer-only arithmetic. No floats.
//!
//! # Respiratory Cycle
//!
//! ```text
//! Inspiration ──(volume reaches tidal target)──→ EndInspiratory
//!      ↑                                              │
//!      │                                   (pause expires)
//!      │                                              ↓
//! EndExpiratory ←──(volume reaches FRC)── Expiration
//!//!  (CPG drive reaches threshold)
//!      └──→ Inspiration
//! ```

use std::time::Instant;

/// Phase of the respiratory cycle.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub enum RespiratoryPhase {
    /// Active diaphragm contraction. Lung volume increasing.
    /// Intrathoracic pressure drops, air flows in.
    /// Vagal tone decreases (RSA: HR increases during inspiration).
    Inspiration,
    /// Brief pause at peak volume. Hering-Breuer reflex onset.
    /// Lung stretch receptors maximally active.
    EndInspiratory,
    /// Passive elastic recoil. Lung volume decreasing.
    /// Intrathoracic pressure rises, air flows out.
    /// Vagal tone returns (RSA: HR decreases during expiration).
    Expiration,
    /// Pause at functional residual capacity (FRC). Quiet before next breath.
    /// CO2 continues to accumulate from metabolism.
    /// CPG drive ramp accumulates toward next inspiration threshold.
    EndExpiratory,
}

/// Lung volume and airflow state.
///
/// All volumes in "lung units" (0-255), where 0 = residual volume
/// and 255 = total lung capacity.
///
/// FRC (functional residual capacity) ~75: the resting volume between breaths.
/// Normal tidal volume ~30 units: breathing oscillates between ~75 and ~105.
pub struct LungState {
    /// Current lung volume (0-255 lung units).
    pub volume: u8,
    /// Current instantaneous airflow. Positive = inspiratory, negative = expiratory.
    pub flow: i16,
    /// Current respiratory phase.
    pub phase: RespiratoryPhase,
    /// When the current phase started.
    pub phase_start: Instant,
    /// Effective tidal volume for this breath cycle (autonomic-modulated).
    pub current_tidal_volume: u8,
}

impl LungState {
    /// Create a new lung state at rest (end-expiratory, FRC volume).
    pub fn new(config: &RespiratoryConfig, now: Instant) -> Self {
        Self {
            volume: config.frc,
            flow: 0,
            phase: RespiratoryPhase::EndExpiratory,
            phase_start: now,
            current_tidal_volume: config.base_tidal_volume,
        }
    }

    /// Compute phase progress as 0-255 (0 = phase just started, 255 = phase ending).
    pub fn phase_progress_256(&self, now: Instant, config: &RespiratoryConfig) -> u8 {
        let elapsed_us = now.duration_since(self.phase_start).as_micros() as u64;
        let phase_duration_us = match self.phase {
            RespiratoryPhase::Inspiration => config.inspiration_us,
            RespiratoryPhase::EndInspiratory => config.end_inspiratory_pause_us,
            RespiratoryPhase::Expiration => config.expiration_us,
            RespiratoryPhase::EndExpiratory => config.end_expiratory_pause_us,
        };
        if phase_duration_us == 0 {
            return 255;
        }
        let progress = (elapsed_us * 255) / phase_duration_us;
        progress.min(255) as u8
    }

    /// Update lung state: advance phase, compute volume and flow.
    ///
    /// Returns `true` if a full breath cycle completed (EndExpiratory → Inspiration transition).
    pub fn update(
        &mut self,
        now: Instant,
        generator: &mut RespiratoryGenerator,
        config: &RespiratoryConfig,
    ) -> bool {
        let elapsed_us = now.duration_since(self.phase_start).as_micros() as u64;

        match self.phase {
            RespiratoryPhase::Inspiration => {
                // Volume rises from FRC toward FRC + tidal volume
                let target_volume = config.frc.saturating_add(self.current_tidal_volume);
                let volume_range = self.current_tidal_volume as u64;
                if config.inspiration_us > 0 && volume_range > 0 {
                    let progress = (elapsed_us * volume_range) / config.inspiration_us;
                    self.volume = config.frc.saturating_add(progress.min(volume_range) as u8);
                    // Flow = volume change rate (units per second)
                    let flow = (volume_range * 1_000_000) / config.inspiration_us;
                    self.flow = flow.min(i16::MAX as u64) as i16;
                }

                if elapsed_us >= config.inspiration_us {
                    self.volume = target_volume;
                    self.flow = 0;
                    self.phase = RespiratoryPhase::EndInspiratory;
                    self.phase_start = now;
                }
                false
            }
            RespiratoryPhase::EndInspiratory => {
                // Brief pause at peak volume
                self.flow = 0;
                if elapsed_us >= config.end_inspiratory_pause_us {
                    self.phase = RespiratoryPhase::Expiration;
                    self.phase_start = now;
                }
                false
            }
            RespiratoryPhase::Expiration => {
                // Volume falls from FRC + tidal back to FRC
                let peak_volume = config.frc.saturating_add(self.current_tidal_volume);
                let volume_range = self.current_tidal_volume as u64;
                if config.expiration_us > 0 && volume_range > 0 {
                    let progress = (elapsed_us * volume_range) / config.expiration_us;
                    let descent = progress.min(volume_range) as u8;
                    self.volume = peak_volume.saturating_sub(descent);
                    // Negative flow (expiratory)
                    let flow = (volume_range * 1_000_000) / config.expiration_us;
                    self.flow = -(flow.min(i16::MAX as u64) as i16);
                }

                if elapsed_us >= config.expiration_us {
                    self.volume = config.frc;
                    self.flow = 0;
                    self.phase = RespiratoryPhase::EndExpiratory;
                    self.phase_start = now;
                    generator.reset_drive();
                }
                false
            }
            RespiratoryPhase::EndExpiratory => {
                // CPG drive ramps. When threshold reached, start next inspiration.
                self.flow = 0;
                self.volume = config.frc;

                // Compute drive from elapsed time (not accumulated)
                let drive = (generator.drive_rate_per_sec as u64 * elapsed_us) / 1_000_000;
                generator.drive = drive.min(255) as u8;

                // Also wait for minimum pause duration
                if generator.drive >= generator.effective_threshold
                    && elapsed_us >= config.end_expiratory_pause_us
                {
                    self.phase = RespiratoryPhase::Inspiration;
                    self.phase_start = now;
                    generator.reset_drive();
                    return true; // New breath cycle starting
                }
                false
            }
        }
    }
}

/// Central pattern generator for respiratory rhythm.
///
/// The pre-Botzinger complex fires rhythmically. Between breaths,
/// a drive ramp accumulates (like the SA node's diastolic depolarization).
/// When drive reaches threshold, inspiration begins.
pub struct RespiratoryGenerator {
    /// Drive accumulator (0-255). Ramps during EndExpiratory.
    pub drive: u8,
    /// Current effective drive rate per second (after autonomic modulation).
    pub drive_rate_per_sec: u32,
    /// Current effective threshold for inspiration (CO2/O2 modulated).
    pub effective_threshold: u8,
    /// Base drive rate (from config, for modulation reference).
    base_drive_rate_per_sec: u32,
}

impl RespiratoryGenerator {
    /// Create a new generator from config.
    pub fn new(config: &RespiratoryConfig) -> Self {
        Self {
            drive: 0,
            drive_rate_per_sec: config.base_drive_rate_per_sec,
            effective_threshold: config.base_drive_threshold,
            base_drive_rate_per_sec: config.base_drive_rate_per_sec,
        }
    }

    /// Reset drive accumulator at the start of a new breath cycle.
    pub fn reset_drive(&mut self) {
        self.drive = 0;
    }

    /// Apply autonomic and chemoreceptor modulation.
    ///
    /// NE increases drive rate (faster breathing). ACh decreases it.
    /// CO2 above baseline lowers threshold (breathe sooner).
    /// O2 below hypoxic threshold provides additional drive.
    pub fn apply_modulation(
        &mut self,
        ne: u8,
        ach: u8,
        co2: u8,
        o2: u8,
        config: &RespiratoryConfig,
    ) {
        let base = self.base_drive_rate_per_sec as u64;

        // NE: up to 2x drive rate at NE=255 with full sensitivity
        let ne_boost = (base * ne as u64 * config.ne_sensitivity as u64) / (255 * 255);

        // ACh: down to 0.5x drive rate at ACh=255 with full sensitivity
        let ach_brake = (base * ach as u64 * config.ach_sensitivity as u64) / (255 * 255 * 2);

        let modulated = (base as i64) + (ne_boost as i64) - (ach_brake as i64);
        let floor = (base / 4).max(1) as i64; // lungs cannot stop
        let ceiling = (base * 4) as i64;
        self.drive_rate_per_sec = modulated.clamp(floor, ceiling) as u32;

        // CO2 chemoreceptor: lower threshold when CO2 above baseline
        let co2_excess = co2.saturating_sub(config.co2_baseline);
        let co2_drive = (co2_excess as u32 * config.co2_sensitivity as u32) / 255;

        // O2 hypoxic drive: additional threshold reduction when O2 drops below threshold
        let o2_drive = if o2 < config.o2_hypoxic_threshold {
            let deficit = config.o2_hypoxic_threshold.saturating_sub(o2);
            (deficit as u32 * config.o2_sensitivity as u32) / 255
        } else {
            0
        };

        // Effective threshold = base - CO2 drive - O2 drive
        let threshold = config.base_drive_threshold as i32
            - co2_drive as i32
            - o2_drive as i32;
        self.effective_threshold = threshold
            .clamp(config.min_drive_threshold as i32, config.base_drive_threshold as i32)
            as u8;
    }
}

/// Gas exchange environment for the lungs.
///
/// CO2 and O2 are tracked as u8 concentrations (0-255).
/// CO2 accumulates continuously from metabolism. O2 depletes continuously.
/// Ventilation (during Inspiration/Expiration phases) clears CO2 and
/// replenishes O2. Clearance scales with tidal volume.
pub struct GasPool {
    /// Arterial CO2 analog (0-255). Primary respiratory drive signal.
    pub co2: u8,
    /// Arterial O2 analog (0-255). Secondary drive (hypoxic response).
    pub o2: u8,
    /// Last time gas metabolism was applied.
    last_metabolism: Instant,
}

impl GasPool {
    /// Create a new gas pool at baseline values.
    pub fn new(config: &RespiratoryConfig, now: Instant) -> Self {
        Self {
            co2: config.co2_baseline,
            o2: config.o2_baseline,
            last_metabolism: now,
        }
    }

    /// Metabolize gases and optionally ventilate.
    ///
    /// CO2 production and O2 consumption happen continuously (metabolism never stops).
    /// When `ventilating` is true (Inspiration/Expiration phases), CO2 clearance
    /// and O2 replenishment also occur, scaled by tidal volume.
    ///
    /// Uses 250ms batching like cardiac metabolism.
    pub fn metabolize(
        &mut self,
        now: Instant,
        ventilating: bool,
        tidal_volume: u8,
        config: &RespiratoryConfig,
    ) {
        let elapsed_us = now.duration_since(self.last_metabolism).as_micros() as u64;
        if elapsed_us < 250_000 {
            return;
        }
        self.last_metabolism = now;

        // CO2 production (continuous — metabolism never stops)
        let co2_produced =
            (config.co2_production_rate_per_sec as u64 * elapsed_us) / 1_000_000;
        self.co2 = self.co2.saturating_add(co2_produced.min(255) as u8);

        // O2 consumption (continuous)
        let o2_consumed =
            (config.o2_consumption_rate_per_sec as u64 * elapsed_us) / 1_000_000;
        self.o2 = self.o2.saturating_sub(o2_consumed.min(255) as u8);

        // Ventilation: clear CO2 and replenish O2 during active breathing
        if ventilating {
            let tidal_ratio_x256 = if config.base_tidal_volume > 0 {
                (tidal_volume as u64 * 256) / config.base_tidal_volume as u64
            } else {
                256
            };

            // CO2 clearance
            let base_clearance =
                (config.co2_clearance_rate_per_sec as u64 * elapsed_us) / 1_000_000;
            let effective_clearance = (base_clearance * tidal_ratio_x256) / 256;
            if self.co2 > config.co2_baseline {
                let distance = (self.co2 - config.co2_baseline) as u64;
                let cleared = effective_clearance.min(distance);
                self.co2 = self.co2.saturating_sub(cleared as u8);
            }

            // O2 replenishment
            let base_replenishment =
                (config.o2_replenishment_rate_per_sec as u64 * elapsed_us) / 1_000_000;
            let effective_replenishment = (base_replenishment * tidal_ratio_x256) / 256;
            if self.o2 < config.o2_baseline {
                let distance = (config.o2_baseline - self.o2) as u64;
                let replenished = effective_replenishment.min(distance);
                self.o2 = self.o2.saturating_add(replenished as u8);
            }
        }
    }
}

/// Compute tidal volume modulated by autonomic state.
///
/// NE deepens breaths. ACh makes them shallower. CO2 excess deepens breaths.
pub fn compute_tidal_volume(ne: u8, ach: u8, co2: u8, config: &RespiratoryConfig) -> u8 {
    let base = config.base_tidal_volume as i32;

    // NE deepens: up to +base at NE=255 with full sensitivity
    let ne_depth = (base * ne as i32 * config.ne_depth_sensitivity as i32) / (255 * 255);

    // ACh shallows: down to -base/2 at ACh=255 with full sensitivity
    let ach_shallow = (base * ach as i32 * config.ach_depth_sensitivity as i32) / (255 * 255 * 2);

    // CO2 deepens: up to +base/2 at max excess
    let co2_excess = co2.saturating_sub(config.co2_baseline) as i32;
    let co2_depth = (base * co2_excess * config.co2_sensitivity as i32) / (255 * 255);

    let modulated = base + ne_depth - ach_shallow + co2_depth;
    modulated.clamp(config.min_tidal_volume as i32, config.max_tidal_volume as i32) as u8
}

/// Complete respiratory system configuration.
///
/// Default values produce resting breathing at ~15 breaths/min with
/// CO2/O2 homeostasis, tidal volume ~30 lung units, and RSA coupling.
#[derive(Clone, Debug)]
pub struct RespiratoryConfig {
    // === Timing ===
    /// Base inspiration duration in microseconds.
    pub inspiration_us: u64,
    /// End-inspiratory pause duration (Hering-Breuer).
    pub end_inspiratory_pause_us: u64,
    /// Base expiration duration in microseconds.
    pub expiration_us: u64,
    /// End-expiratory pause base duration.
    pub end_expiratory_pause_us: u64,

    // === Volume ===
    /// Functional residual capacity — resting lung volume (0-255).
    pub frc: u8,
    /// Resting tidal volume (0-255 lung units).
    pub base_tidal_volume: u8,
    /// Maximum tidal volume (deep breath).
    pub max_tidal_volume: u8,
    /// Minimum tidal volume (shallow breathing).
    pub min_tidal_volume: u8,

    // === CPG Drive ===
    /// Base drive rate per second.
    pub base_drive_rate_per_sec: u32,
    /// Base drive threshold for initiating inspiration.
    pub base_drive_threshold: u8,
    /// Minimum threshold (under high CO2 / low O2).
    pub min_drive_threshold: u8,

    // === Autonomic Sensitivity ===
    /// NE sensitivity for respiratory rate modulation (0-255).
    pub ne_sensitivity: u8,
    /// ACh sensitivity for respiratory rate modulation (0-255).
    pub ach_sensitivity: u8,
    /// NE sensitivity for tidal volume modulation (0-255).
    pub ne_depth_sensitivity: u8,
    /// ACh sensitivity for tidal volume modulation (0-255).
    pub ach_depth_sensitivity: u8,

    // === Chemoreceptor ===
    /// CO2 baseline (homeostatic target).
    pub co2_baseline: u8,
    /// O2 baseline.
    pub o2_baseline: u8,
    /// CO2 sensitivity for threshold reduction (0-255).
    pub co2_sensitivity: u8,
    /// O2 threshold below which hypoxic drive kicks in.
    pub o2_hypoxic_threshold: u8,
    /// O2 sensitivity for threshold reduction (0-255).
    pub o2_sensitivity: u8,

    // === Gas Exchange ===
    /// Resting CO2 production rate (units/sec).
    pub co2_production_rate_per_sec: u32,
    /// Resting O2 consumption rate (units/sec).
    pub o2_consumption_rate_per_sec: u32,
    /// Base CO2 clearance rate during ventilation (units/sec).
    pub co2_clearance_rate_per_sec: u32,
    /// Base O2 replenishment rate during ventilation (units/sec).
    pub o2_replenishment_rate_per_sec: u32,

    // === RSA Coupling ===
    /// Whether RSA coupling to the heart is active.
    pub rsa_enabled: bool,
    /// Resting vagal tone that RSA modulates around (ACh amount).
    pub rsa_vagal_baseline: u8,
    /// Depth of RSA modulation (0-255). How much vagal tone varies.
    pub rsa_depth: u8,
    /// Extra vagal augmentation during late expiration.
    pub rsa_expiratory_augmentation: u8,
}

impl Default for RespiratoryConfig {
    fn default() -> Self {
        Self {
            // ~4 second cycle = ~15 breaths/min
            inspiration_us: 1_600_000,
            end_inspiratory_pause_us: 100_000,
            expiration_us: 2_200_000,
            end_expiratory_pause_us: 100_000,

            frc: 75,
            base_tidal_volume: 30,
            max_tidal_volume: 120,
            min_tidal_volume: 10,

            base_drive_rate_per_sec: 64,
            base_drive_threshold: 128,
            min_drive_threshold: 40,

            ne_sensitivity: 128,
            ach_sensitivity: 96,
            ne_depth_sensitivity: 100,
            ach_depth_sensitivity: 80,

            co2_baseline: 40,
            o2_baseline: 200,
            co2_sensitivity: 128,
            o2_hypoxic_threshold: 120,
            o2_sensitivity: 64,

            co2_production_rate_per_sec: 20,
            o2_consumption_rate_per_sec: 15,
            // Clearance must compensate for production during non-ventilating
            // phases (EndInspiratory + EndExpiratory ≈ 36% of cycle).
            // Need: clearance ≥ production / ventilation_fraction ≈ 20/0.64 ≈ 31.
            // Use 35 for comfortable margin.
            co2_clearance_rate_per_sec: 35,
            o2_replenishment_rate_per_sec: 28,

            rsa_enabled: true,
            rsa_vagal_baseline: 30,
            rsa_depth: 180,
            rsa_expiratory_augmentation: 15,
        }
    }
}

/// Compute RSA vagal modulation based on breath phase.
///
/// Returns the ACh amount (0-255) to send to the heart via the RSA signal.
/// Replaces the heart's ACh baseline with a breath-phase-modulated signal.
///
/// Inspiration: vagal tone decreases (HR rises).
/// Expiration: vagal tone increases (HR falls).
pub fn compute_vagal_modulation(
    phase: RespiratoryPhase,
    phase_progress_256: u8,
    config: &RespiratoryConfig,
) -> u8 {
    let baseline = config.rsa_vagal_baseline;
    let depth = config.rsa_depth;

    // Maximum vagal reduction during peak inspiration
    let max_reduction = (baseline as u16 * depth as u16) / 255;
    let min_vagal = baseline.saturating_sub(max_reduction as u8);
    let augmented = baseline.saturating_add(config.rsa_expiratory_augmentation);

    match phase {
        RespiratoryPhase::Inspiration => {
            // Vagal withdrawal: ACh decreases as inspiration progresses.
            // Linear ramp from baseline to min_vagal.
            let reduction =
                (max_reduction * phase_progress_256 as u16) / 255;
            baseline.saturating_sub(reduction as u8)
        }
        RespiratoryPhase::EndInspiratory => {
            // Minimal vagal tone (peak stretch).
            min_vagal
        }
        RespiratoryPhase::Expiration => {
            // Vagal return: linear ramp from min_vagal to augmented.
            let range = augmented.saturating_sub(min_vagal) as u16;
            let current =
                min_vagal as u16 + (range * phase_progress_256 as u16) / 255;
            current.min(255) as u8
        }
        RespiratoryPhase::EndExpiratory => {
            // Full vagal tone restored + augmentation.
            augmented
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use std::time::Duration;

    #[test]
    fn phase_transitions_are_sequential() {
        let config = RespiratoryConfig::default();
        let base = Instant::now();
        let mut lung = LungState::new(&config, base);
        let mut gen = RespiratoryGenerator::new(&config);

        // Start in EndExpiratory. Force CPG to fire immediately.
        gen.effective_threshold = 1;
        gen.drive_rate_per_sec = 10000; // Very fast ramp

        let mut phases = vec![lung.phase];

        // Advance through at least one full cycle
        for ms in 1..10000u64 {
            let now = base + Duration::from_millis(ms);
            lung.update(now, &mut gen, &config);
            if *phases.last().unwrap() != lung.phase {
                phases.push(lung.phase);
            }
            if phases.len() >= 6 {
                break;
            }
        }

        // Should see: EndExpiratory → Inspiration → EndInspiratory → Expiration → EndExpiratory → Inspiration
        assert!(phases.len() >= 5, "Got phases: {:?}", phases);
        assert_eq!(phases[0], RespiratoryPhase::EndExpiratory);
        assert_eq!(phases[1], RespiratoryPhase::Inspiration);
        assert_eq!(phases[2], RespiratoryPhase::EndInspiratory);
        assert_eq!(phases[3], RespiratoryPhase::Expiration);
        assert_eq!(phases[4], RespiratoryPhase::EndExpiratory);
    }

    #[test]
    fn volume_rises_during_inspiration() {
        let config = RespiratoryConfig::default();
        let base = Instant::now();
        let mut lung = LungState::new(&config, base);
        let mut gen = RespiratoryGenerator::new(&config);

        // Force into Inspiration
        lung.phase = RespiratoryPhase::Inspiration;
        lung.phase_start = base;
        lung.volume = config.frc;

        // Advance 800ms into a 1600ms inspiration
        let now = base + Duration::from_millis(800);
        lung.update(now, &mut gen, &config);

        assert!(lung.volume > config.frc, "Volume should rise during inspiration");
        assert!(lung.flow > 0, "Flow should be positive (inspiratory)");
    }

    #[test]
    fn volume_falls_during_expiration() {
        let config = RespiratoryConfig::default();
        let base = Instant::now();
        let mut lung = LungState::new(&config, base);
        let mut gen = RespiratoryGenerator::new(&config);

        // Force into Expiration at peak volume
        let peak = config.frc.saturating_add(config.base_tidal_volume);
        lung.phase = RespiratoryPhase::Expiration;
        lung.phase_start = base;
        lung.volume = peak;

        // Advance 1100ms into a 2200ms expiration
        let now = base + Duration::from_millis(1100);
        lung.update(now, &mut gen, &config);

        assert!(lung.volume < peak, "Volume should fall during expiration");
        assert!(lung.flow < 0, "Flow should be negative (expiratory)");
    }

    #[test]
    fn cpg_drive_accumulates_during_pause() {
        let config = RespiratoryConfig::default();
        let base = Instant::now();
        let mut lung = LungState::new(&config, base);
        let mut gen = RespiratoryGenerator::new(&config);

        // In EndExpiratory, drive should ramp up
        assert_eq!(lung.phase, RespiratoryPhase::EndExpiratory);
        assert_eq!(gen.drive, 0);

        // Advance 1 second (drive_rate=64 → drive=64)
        let now = base + Duration::from_secs(1);
        lung.update(now, &mut gen, &config);

        assert!(gen.drive > 0, "CPG drive should accumulate during EndExpiratory");
        // At 64/sec for 1 second = 64
        assert_eq!(gen.drive, 64);
    }

    #[test]
    fn co2_lowers_cpg_threshold() {
        let config = RespiratoryConfig::default();
        let mut gen = RespiratoryGenerator::new(&config);

        // At baseline CO2: threshold should be base
        gen.apply_modulation(0, 0, config.co2_baseline, config.o2_baseline, &config);
        assert_eq!(gen.effective_threshold, config.base_drive_threshold);

        // At elevated CO2: threshold should drop
        gen.apply_modulation(0, 0, config.co2_baseline + 50, config.o2_baseline, &config);
        assert!(
            gen.effective_threshold < config.base_drive_threshold,
            "High CO2 should lower threshold: got {}",
            gen.effective_threshold
        );
    }

    #[test]
    fn ne_increases_drive_rate() {
        let config = RespiratoryConfig::default();
        let mut gen = RespiratoryGenerator::new(&config);

        gen.apply_modulation(0, 0, config.co2_baseline, config.o2_baseline, &config);
        let resting_rate = gen.drive_rate_per_sec;

        gen.apply_modulation(200, 0, config.co2_baseline, config.o2_baseline, &config);
        assert!(
            gen.drive_rate_per_sec > resting_rate,
            "NE should increase drive rate"
        );
    }

    #[test]
    fn ach_decreases_drive_rate() {
        let config = RespiratoryConfig::default();
        let mut gen = RespiratoryGenerator::new(&config);

        gen.apply_modulation(0, 0, config.co2_baseline, config.o2_baseline, &config);
        let resting_rate = gen.drive_rate_per_sec;

        gen.apply_modulation(0, 200, config.co2_baseline, config.o2_baseline, &config);
        assert!(
            gen.drive_rate_per_sec < resting_rate,
            "ACh should decrease drive rate"
        );
    }

    #[test]
    fn gas_exchange_clears_co2_during_ventilation() {
        let config = RespiratoryConfig::default();
        let base = Instant::now();
        let mut gas = GasPool::new(&config, base);

        // Elevate CO2 above baseline
        gas.co2 = 80;

        // Run metabolism + ventilation for 500ms
        let now = base + Duration::from_millis(500);
        gas.metabolize(now, true, config.base_tidal_volume, &config);

        // CO2 should have decreased (clearance > production at elevated levels)
        assert!(
            gas.co2 < 80,
            "Ventilation should clear CO2: got {}",
            gas.co2
        );
    }

    #[test]
    fn co2_accumulates_during_pause() {
        let config = RespiratoryConfig::default();
        let base = Instant::now();
        let mut gas = GasPool::new(&config, base);

        // Only run metabolism (no ventilation) for 500ms
        let now = base + Duration::from_millis(500);
        gas.metabolize(now, false, 0, &config);

        // CO2 should have risen above baseline
        assert!(
            gas.co2 > config.co2_baseline,
            "CO2 should rise without ventilation: got {} vs baseline {}",
            gas.co2,
            config.co2_baseline
        );
    }
}