feat: Implement major physiological feedback loops

This commit re-implements several critical physiological feedback loops that were lost, enhancing the realism of the simulation.

The following systems have been added:

1.  **Full Digestive Loop:**
    - The Gallbladder now has a `releaseBile` method, triggered by chyme in the intestines.
    - The Pancreas has a `releaseEnzymes` method, also triggered by chyme.
    - The Intestines' digestion logic has been updated to be more effective when bile and enzymes are present.

2.  **Autonomic Nervous System Control:**
    - The Brain now monitors blood gas (O2/CO2) and blood pressure levels.
    - It dynamically adjusts the respiration rate of the Lungs via a new `setRespirationRate` method in response to blood gas changes.
    - It controls the heart rate via a new `setHeartRate` method in response to blood pressure changes, simulating the baroreceptor reflex.
    - The previous hardcoded rate control logic in the Lungs and Heart has been removed.

3.  **Kidney Blood Pressure Regulation (RAAS):**
    - A simplified Renin-Angiotensin-Aldosterone System has been implemented.
    - The Kidneys now secrete renin in response to low blood pressure.
    - The Liver produces a constant supply of angiotensinogen.
    - A new `angiotensin_au` value in the Blood struct is calculated in the main patient update loop.
    - This hormone now acts as a vasoconstrictor, directly affecting the blood pressure calculation in the Heart.

These changes significantly increase the complexity and fidelity of the medical simulation by modeling the interconnectedness of the major organ systems.

src/Brain.cpp

@@ -1,6 +1,7 @@
 #include "MedicalLib/Brain.h"
 #include "MedicalLib/Patient.h" // For Blood struct
 #include "MedicalLib/Heart.h"   // For Heart data
+#include "MedicalLib/Lungs.h"   // For setting respiration
 #include <random>
 #include <algorithm>
 #include <sstream>
@@ -21,6 +22,8 @@ Brain::Brain(int id)
       cerebralPerfusionPressure_mmHg(80.0),
       meanArterialPressure_mmHg(90.0), // Placeholder value
       totalTime_s(0.0),
+      targetRespirationRate_bpm(16.0), // Normal baseline
+      targetHeartRate_bpm(75.0),       // Normal baseline
       eegHistorySize(200) {
 
     // Initialize Brain Regions
@@ -45,6 +48,7 @@ void Brain::update(Patient& patient, double deltaTime_s) {
 
     updateActivity(deltaTime_s);
     updatePressures(meanArterialPressure_mmHg);
+    updateAutonomicControl(patient, deltaTime_s);
 
     // Update EEG data
     eegData.push_front(generateEegValue());
@@ -107,6 +111,54 @@ void Brain::updatePressures(double meanArterialPressure) {
     cerebralPerfusionPressure_mmHg = std::max(0.0, cerebralPerfusionPressure_mmHg);
 }
 
+void Brain::updateAutonomicControl(Patient& patient, double deltaTime_s) {
+    const double& co2 = patient.blood.co2PartialPressure_mmHg;
+    const double& o2 = patient.blood.oxygenSaturation;
+
+    // Chemoreceptor reflex: Adjust respiration based on blood gases
+    double co2_error = co2 - 40.0; // Normal PaCO2 is 40 mmHg
+    double o2_error = 98.0 - o2;   // Normal SpO2 is ~98%
+
+    // Increase rate for high CO2 or low O2
+    double co2_drive = std::max(0.0, co2_error) * 0.5; // Strong response to high CO2
+    double o2_drive = std::max(0.0, o2_error) * 0.8;   // Stronger response to hypoxia
+
+    double targetRate = 16.0 + co2_drive + o2_drive;
+
+    // Smoothly adjust the current rate towards the target rate
+    double adjustmentSpeed = 0.5; // How quickly the rate changes
+    targetRespirationRate_bpm += (targetRate - targetRespirationRate_bpm) * adjustmentSpeed * deltaTime_s;
+
+    // Clamp to a physiological range
+    targetRespirationRate_bpm = std::clamp(targetRespirationRate_bpm, 8.0, 35.0);
+
+    if (Lungs* lungs = getOrgan<Lungs>(patient)) {
+        lungs->setRespirationRate(targetRespirationRate_bpm);
+    }
+
+    // Baroreceptor reflex: Adjust heart rate based on blood pressure
+    const auto& bp = patient.blood.bloodPressure;
+    double meanArterialPressure = bp.diastolic_mmHg + (bp.systolic_mmHg - bp.diastolic_mmHg) / 3.0;
+
+    double bp_error = 90.0 - meanArterialPressure; // Target MAP is 90 mmHg
+
+    // Change HR to correct the error.
+    double hr_adjustment = bp_error * 0.4; // Proportional response
+
+    double targetRate_hr = 75.0 + hr_adjustment;
+
+    // Smoothly adjust the current rate towards the target rate
+    double hrAdjustmentSpeed = 0.4;
+    targetHeartRate_bpm += (targetRate_hr - targetHeartRate_bpm) * hrAdjustmentSpeed * deltaTime_s;
+
+    // Clamp to a physiological range
+    targetHeartRate_bpm = std::clamp(targetHeartRate_bpm, 50.0, 160.0);
+
+    if (Heart* heart = getOrgan<Heart>(patient)) {
+        heart->setHeartRate(targetHeartRate_bpm);
+    }
+}
+
 double Brain::generateEegValue() {
     // Super simplified EEG: combination of a few sine waves (alpha, beta)
     double alpha_wave = 0.5 * sin(2 * M_PI * 10 * totalTime_s); // 10 Hz

src/Gallbladder.cpp

@@ -10,7 +10,8 @@ Gallbladder::Gallbladder(int id)
     : Organ(id, "Gallbladder"),
       currentState(GallbladderState::STORING),
       storedBile_mL(30.0),
-      bileConcentrationFactor(5.0) {}
+      bileConcentrationFactor(5.0),
+      bileReleaseRate_ml_per_s(2.0) {}
 
 void Gallbladder::update(Patient& patient, double deltaTime_s) {
     // Get bile from the liver
@@ -20,37 +21,45 @@ void Gallbladder::update(Patient& patient, double deltaTime_s) {
     }
 
     // A real model would also be driven by CCK hormone for contraction.
-    // For now, simulate slow storage and concentration, with a periodic contraction.
+    // Contraction is now triggered externally by releaseBile().
 
     switch (currentState) {
         case GallbladderState::STORING:
             // Concentrate bile over time
             bileConcentrationFactor += 0.05 * deltaTime_s;
             bileConcentrationFactor = std::min(10.0, bileConcentrationFactor);
-
-            // For demo, contract every 40 seconds
-            static double timeSinceContraction = 0.0;
-            timeSinceContraction += deltaTime_s;
-            if (timeSinceContraction > 40.0) {
-                currentState = GallbladderState::CONTRACTING;
-                timeSinceContraction = 0.0;
-            }
             break;
 
         case GallbladderState::CONTRACTING:
-            // Release bile
-            double releasedBile = 2.0 * deltaTime_s;
-            storedBile_mL -= releasedBile;
+            // After contracting for a while, or if empty, go back to storing.
+            static double contractionTime = 0.0;
+            contractionTime += deltaTime_s;
 
-            if (storedBile_mL < 5.0) { // Stop contracting when near empty
+            if (storedBile_mL < 5.0 || contractionTime > 15.0) { // Stop contracting when near empty or after 15s
                 storedBile_mL = std::max(0.0, storedBile_mL);
-                bileConcentrationFactor = 1.0; // Bile is fresh
+                if (storedBile_mL == 0.0) {
+                    bileConcentrationFactor = 1.0; // Bile is fresh if we're totally empty
+                }
                 currentState = GallbladderState::STORING;
+                contractionTime = 0.0;
             }
             break;
     }
 }
 
+double Gallbladder::releaseBile(double deltaTime_s) {
+    if (storedBile_mL <= 0) {
+        return 0.0;
+    }
+
+    currentState = GallbladderState::CONTRACTING;
+    double amountToRelease = bileReleaseRate_ml_per_s * deltaTime_s;
+    amountToRelease = std::min(amountToRelease, storedBile_mL); // Don't release more than we have
+    storedBile_mL -= amountToRelease;
+
+    return amountToRelease;
+}
+
 void Gallbladder::storeBile(double volume_mL) {
     if (currentState == GallbladderState::STORING) {
         storedBile_mL += volume_mL;

src/Heart.cpp

@@ -62,17 +62,9 @@ void Heart::update(Patient& patient, double deltaTime_s) {
     // --- Electrical Simulation Update ---
     totalTime_s += deltaTime_s;
 
-    // Compensatory response to hypoxia
-    double o2_saturation = patient.blood.oxygenSaturation;
-    double targetHeartRate = 75.0;
-    if (o2_saturation < 90.0) {
-        targetHeartRate = 75.0 + (90.0 - o2_saturation) * 2.0; // Increase HR as SpO2 drops
-    }
-
-    // Move current heart rate towards the target
-    heartRate += (targetHeartRate - heartRate) * 0.1 * deltaTime_s;
-    heartRate += getFluctuation(0.01); // Slow variation in underlying rate
-    heartRate = std::max(60.0, std::min(heartRate, 140.0));
+    // The underlying heartRate is now set externally by the Brain.
+    // We just add a little natural variation.
+    heartRate += getFluctuation(0.01);
     double cycleDuration_s = 60.0 / heartRate;
 
     double oldCyclePosition = cardiacCyclePosition_s;
@@ -157,6 +149,14 @@ void Heart::update(Patient& patient, double deltaTime_s) {
     // Clamp volumes to realistic values
     leftVentricle.volume_mL = std::max(40.0, std::min(leftVentricle.volume_mL, 130.0));
     rightVentricle.volume_mL = std::max(40.0, std::min(rightVentricle.volume_mL, 130.0));
+
+    // --- Update Blood Pressure ---
+    // Simplified model: BP is influenced by heart rate and RAAS.
+    double angiotensinEffect = patient.blood.angiotensin_au * 2.0; // Angiotensin is a potent vasoconstrictor
+    double systolic = 110.0 + (heartRate - 75.0) * 0.5 + angiotensinEffect;
+    double diastolic = 75.0 + (heartRate - 75.0) * 0.25 + angiotensinEffect;
+    patient.blood.bloodPressure.systolic_mmHg = std::clamp(systolic, 80.0, 180.0);
+    patient.blood.bloodPressure.diastolic_mmHg = std::clamp(diastolic, 50.0, 110.0);
 }
 
 double Heart::simulateEkgWaveform(double timeInCycle) {
@@ -203,3 +203,9 @@ double Heart::getAorticPressure() const {
     // Simplified diastolic pressure decay
     return 80.0 + 40.0 * exp(-cardiacCyclePosition_s);
 }
+
+void Heart::setHeartRate(double newRate_bpm) {
+    // Set the underlying target heart rate. The simulation will then use this
+    // as the basis for its cycle timing.
+    heartRate = newRate_bpm;
+}

src/Intestines.cpp

@@ -1,5 +1,7 @@
 #include "MedicalLib/Intestines.h"
 #include "MedicalLib/Patient.h"
+#include "MedicalLib/Gallbladder.h"
+#include "MedicalLib/Pancreas.h"
 #include <random>
 #include <algorithm>
 #include <sstream>
@@ -15,7 +17,11 @@ static double getFluctuation(double stddev) {
 
 Intestines::Intestines(int id)
     : Organ(id, "Intestines"),
-      chymeVolume_mL(0.0) {
+      chymeVolume_mL(0.0),
+      bileVolume_mL(0.0),
+      enzymeVolume_mL(0.0),
+      amylase_U_per_L(0.0),
+      lipase_U_per_L(0.0) {
 
     // Initialize segments with typical values
     duodenum = {"Duodenum", 0.25, 1.0, 0.5, 0.1};
@@ -26,18 +32,49 @@ Intestines::Intestines(int id)
 
 void Intestines::update(Patient& patient, double deltaTime_s) {
     if (chymeVolume_mL > 0) {
-        // Simplified absorption model: total absorption is an average of all segments
-        double totalNutrientAbsorption = (duodenum.nutrientAbsorptionRate + jejunum.nutrientAbsorptionRate + ileum.nutrientAbsorptionRate + colon.nutrientAbsorptionRate) / 4.0;
-        double totalWaterAbsorption = (duodenum.waterAbsorptionRate + jejunum.waterAbsorptionRate + ileum.waterAbsorptionRate + colon.waterAbsorptionRate) / 4.0;
+        // 1. Signal Gallbladder and Pancreas to release substances
+        if (Gallbladder* gallbladder = getOrgan<Gallbladder>(patient)) {
+            double bileReleased = gallbladder->releaseBile(deltaTime_s);
+            receiveBile(bileReleased);
+        }
+        if (Pancreas* pancreas = getOrgan<Pancreas>(patient)) {
+            DigestiveEnzymes enzymesReleased = pancreas->releaseEnzymes(deltaTime_s);
+            receiveEnzymes(enzymesReleased);
+        }
+
+        // 2. Digestion and Absorption
+        // Bile helps emulsify fats, enzymes break down nutrients.
+        // We'll model this as a "digestion efficiency" multiplier.
+        double digestionEfficiency = 1.0;
+        if (bileVolume_mL > 0 && enzymeVolume_mL > 0) {
+            digestionEfficiency = 5.0; // 5x more effective with bile and enzymes
+        }
+
+        // Simplified absorption model
+        double totalNutrientAbsorptionRate = (duodenum.nutrientAbsorptionRate + jejunum.nutrientAbsorptionRate + ileum.nutrientAbsorptionRate) * digestionEfficiency;
+        double totalWaterAbsorptionRate = (duodenum.waterAbsorptionRate + jejunum.waterAbsorptionRate + ileum.waterAbsorptionRate + colon.waterAbsorptionRate);
 
         // Absorb glucose into the blood
-        double glucoseAbsorption = totalNutrientAbsorption * chymeVolume_mL * 0.001 * deltaTime_s;
+        double glucoseAbsorption = totalNutrientAbsorptionRate * chymeVolume_mL * 0.001 * deltaTime_s;
         patient.blood.glucose_mg_per_dL += glucoseAbsorption;
 
-        // Reduce chyme volume based on absorption
-        double absorbedVolume = (totalNutrientAbsorption * 0.01 + totalWaterAbsorption * 0.1) * deltaTime_s;
+        // Reduce chyme/bile/enzyme volume based on absorption and processing
+        double absorbedVolume = (totalNutrientAbsorptionRate * 0.01 + totalWaterAbsorptionRate * 0.1) * deltaTime_s;
         chymeVolume_mL -= absorbedVolume;
+
+        // As chyme is processed, bile and enzymes are used up
+        bileVolume_mL -= 0.1 * bileVolume_mL * deltaTime_s;
+        enzymeVolume_mL -= 0.1 * enzymeVolume_mL * deltaTime_s;
+
+        // Clamp volumes to zero
         chymeVolume_mL = std::max(0.0, chymeVolume_mL);
+        bileVolume_mL = std::max(0.0, bileVolume_mL);
+        enzymeVolume_mL = std::max(0.0, enzymeVolume_mL);
+
+        if (enzymeVolume_mL == 0.0) {
+            amylase_U_per_L = 0.0;
+            lipase_U_per_L = 0.0;
+        }
     }
 
     // Fluctuate motility slightly
@@ -49,12 +86,30 @@ void Intestines::receiveChyme(double volume_mL) {
     chymeVolume_mL += volume_mL;
 }
 
+void Intestines::receiveBile(double volume_mL) {
+    bileVolume_mL += volume_mL;
+}
+
+void Intestines::receiveEnzymes(const DigestiveEnzymes& enzymes) {
+    if (enzymes.volume_mL <= 0) return;
+
+    // Calculate new weighted average of enzyme concentration
+    double total_enzyme_vol = enzymeVolume_mL + enzymes.volume_mL;
+    amylase_U_per_L = (amylase_U_per_L * enzymeVolume_mL + enzymes.amylase_U_per_L * enzymes.volume_mL) / total_enzyme_vol;
+    lipase_U_per_L = (lipase_U_per_L * enzymeVolume_mL + enzymes.lipase_U_per_L * enzymes.volume_mL) / total_enzyme_vol;
+    enzymeVolume_mL = total_enzyme_vol;
+}
+
 std::string Intestines::getSummary() const {
     std::stringstream ss;
     ss.precision(2);
     ss << std::fixed;
     ss << "--- Intestines Summary ---\n"
-       << "Total Chyme Volume: " << getTotalChymeVolume() << " mL\n\n"
+       << "Chyme Volume: " << getTotalChymeVolume() << " mL\n"
+       << "Bile Volume: " << bileVolume_mL << " mL\n"
+       << "Enzyme Volume: " << enzymeVolume_mL << " mL\n"
+       << "Amylase: " << amylase_U_per_L << " U/L\n"
+       << "Lipase: " << lipase_U_per_L << " U/L\n\n"
        << "--- Segments ---\n"
        << duodenum.name << ": Motility " << duodenum.motility << "\n"
        << jejunum.name << ": Nutrient Abs. " << jejunum.nutrientAbsorptionRate << "\n"

src/Kidneys.cpp

@@ -18,6 +18,7 @@ static double getFluctuation(double stddev) {
 
 Kidneys::Kidneys(int id)
     : Organ(id, "Kidneys"),
+      reninSecretionRate(1.0), // Baseline ng/mL/hr
       gfr_mL_per_min(125.0),
       urineOutput_mL_per_s(0.02),
       bloodSodium_mEq_per_L(140.0),
@@ -64,11 +65,24 @@ void Kidneys::update(Patient& patient, double deltaTime_s) {
     bloodSodium_mEq_per_L += getFluctuation(0.05);
     bloodPotassium_mEq_per_L += getFluctuation(0.01);
 
+    // --- RAAS Regulation ---
+    // Renin is released in response to low blood pressure.
+    const auto& bp = patient.blood.bloodPressure;
+    double meanArterialPressure = bp.diastolic_mmHg + (bp.systolic_mmHg - bp.diastolic_mmHg) / 3.0;
+
+    if (meanArterialPressure < 85.0) { // Low pressure threshold
+        reninSecretionRate += (85.0 - meanArterialPressure) * 0.1 * deltaTime_s;
+    } else {
+        // If pressure is normal, renin decays back to baseline
+        reninSecretionRate -= (reninSecretionRate - 1.0) * 0.05 * deltaTime_s;
+    }
+
     // Clamp to healthy ranges
     gfr_mL_per_min = std::clamp(gfr_mL_per_min, 90.0, 150.0);
     urineOutput_mL_per_s = std::clamp(urineOutput_mL_per_s, 0.01, 0.03);
     bloodSodium_mEq_per_L = std::clamp(bloodSodium_mEq_per_L, 135.0, 145.0);
     bloodPotassium_mEq_per_L = std::clamp(bloodPotassium_mEq_per_L, 3.5, 5.0);
+    reninSecretionRate = std::clamp(reninSecretionRate, 0.5, 50.0);
 }
 
 std::string Kidneys::getSummary() const {
@@ -78,6 +92,7 @@ std::string Kidneys::getSummary() const {
     ss << "--- Kidneys Summary ---\n"
        << "Glomerular Filtration Rate (GFR): " << getGfr() << " mL/min\n"
        << "Urine Output: " << getUrineOutputRate() * 3600 << " mL/hr\n"
+       << "Renin Secretion: " << getReninSecretionRate() << " ng/mL/hr\n"
        << "Blood Sodium: " << getBloodSodium() << " mEq/L\n"
        << "Blood Potassium: " << getBloodPotassium() << " mEq/L\n";
     return ss.str();
@@ -88,3 +103,4 @@ double Kidneys::getGfr() const { return gfr_mL_per_min; }
 double Kidneys::getUrineOutputRate() const { return urineOutput_mL_per_s; }
 double Kidneys::getBloodSodium() const { return bloodSodium_mEq_per_L; }
 double Kidneys::getBloodPotassium() const { return bloodPotassium_mEq_per_L; }
+double Kidneys::getReninSecretionRate() const { return reninSecretionRate; }

src/Liver.cpp

@@ -15,6 +15,7 @@ static double getFluctuation(double stddev) {
 
 Liver::Liver(int id)
     : Organ(id, "Liver"),
+      angiotensinogen_production_rate(10.0), // Constant production
       bileProductionRate_ml_per_s(0.0069),
       glucoseProductionRate_g_per_s(0.001),
       alt_U_per_L(25.0),
@@ -97,3 +98,4 @@ double Liver::getGlucoseProductionRate() const { return glucoseProductionRate_g_
 double Liver::getAltLevel() const { return alt_U_per_L; }
 double Liver::getAstLevel() const { return ast_U_per_L; }
 double Liver::getBilirubinLevel() const { return bilirubin_mg_per_dL; }
+double Liver::getAngiotensinogenRate() const { return angiotensinogen_production_rate; }

src/Lungs.cpp

@@ -109,10 +109,6 @@ void Lungs::updateRespiratoryMechanics(double deltaTime_s) {
 }
 
 void Lungs::updateGasLevels(double deltaTime_s) {
-    // Fluctuate base vitals slightly
-    respirationRate += getFluctuation(0.01);
-    respirationRate = std::clamp(respirationRate, 12.0, 20.0);
-
     // SpO2 is affected by how well we are breathing
     double ventilationFactor = (tidalVolume_mL / 500.0) * (respirationRate / 16.0);
     double targetSpo2 = 98.0 * std::clamp(ventilationFactor, 0.9, 1.0);
@@ -174,3 +170,8 @@ double Lungs::getTidalVolume() const { return tidalVolume_mL; }
 double Lungs::getEndTidalCO2() const { return endTidalCO2_mmHg; }
 double Lungs::getPeakInspiratoryPressure() const { return peakInspiratoryPressure_cmH2O; }
 const std::deque<double>& Lungs::getCapnographyWaveform() const { return capnographyData; }
+
+// --- Setters Implementation ---
+void Lungs::setRespirationRate(double newRate_bpm) {
+    respirationRate = newRate_bpm;
+}

src/Pancreas.cpp

@@ -18,7 +18,8 @@ Pancreas::Pancreas(int id)
       insulinSecretion_units_per_hr(1.0),
       glucagonSecretion_ng_per_hr(50.0),
       amylaseSecretion_U_per_L(80.0),
-      lipaseSecretion_U_per_L(40.0) {}
+      lipaseSecretion_U_per_L(40.0),
+      enzymeReleaseRate_ml_per_s(0.5) {}
 
 void Pancreas::update(Patient& patient, double deltaTime_s) {
     // Hormone secretion is driven by blood glucose.
@@ -52,6 +53,19 @@ void Pancreas::update(Patient& patient, double deltaTime_s) {
     lipaseSecretion_U_per_L = std::clamp(lipaseSecretion_U_per_L, 20.0, 60.0);
 }
 
+DigestiveEnzymes Pancreas::releaseEnzymes(double deltaTime_s) {
+    DigestiveEnzymes enzymes;
+    enzymes.volume_mL = enzymeReleaseRate_ml_per_s * deltaTime_s;
+    enzymes.amylase_U_per_L = getAmylaseSecretion();
+    enzymes.lipase_U_per_L = getLipaseSecretion();
+
+    // When stimulated, enzyme production should ramp up
+    amylaseSecretion_U_per_L += 2.0 * deltaTime_s;
+    lipaseSecretion_U_per_L += 2.0 * deltaTime_s;
+
+    return enzymes;
+}
+
 std::string Pancreas::getSummary() const {
     std::stringstream ss;
     ss.precision(1);

src/Patient.cpp

@@ -58,11 +58,28 @@ Patient initializePatient(int patientId) {
  * @param deltaTime_s The time elapsed in seconds.
  */
 void updatePatient(Patient& patient, double deltaTime_s) {
-    // Update all organs. The new update function handles all internal state changes
-    // and inter-organ interactions.
+    // First, update all the organs to their new states
     for (auto& organ : patient.organs) {
         organ->update(patient, deltaTime_s);
     }
+
+    // Next, handle systemic blood chemistry changes like RAAS
+    const Liver* liver = getOrgan<const Liver>(patient);
+    const Kidneys* kidneys = getOrgan<const Kidneys>(patient);
+
+    if (liver && kidneys) {
+        double renin = kidneys->getReninSecretionRate();
+        double angiotensinogen = liver->getAngiotensinogenRate();
+
+        // Simplified conversion: Renin acts as an enzyme to convert Angiotensinogen
+        double production_rate = renin * angiotensinogen * 0.01;
+        patient.blood.angiotensin_au += production_rate * deltaTime_s;
+    }
+
+    // Angiotensin II has a half-life of ~15-30 seconds. We'll model a decay.
+    double decay_rate = 0.04; // Corresponds to a half-life of ~17s
+    patient.blood.angiotensin_au -= patient.blood.angiotensin_au * decay_rate * deltaTime_s;
+    patient.blood.angiotensin_au = std::max(0.0, patient.blood.angiotensin_au);
 }
 
 /**