feat: Implement detailed mechanical heart simulation
This commit introduces a major enhancement to the heart simulation, adding a detailed mechanical model of the cardiac cycle that is synchronized with the existing EKG waveform. Key features: - The Heart class is refactored to include data structures for the four chambers (atria and ventricles) and four valves (mitral, tricuspid, aortic, pulmonary). - The `update` method now simulates the phases of the cardiac cycle, including atrial and ventricular systole and diastole. - Chamber pressures and volumes are dynamically updated based on the phase of the cycle. - Valve statuses (open/closed) are determined by pressure gradients between chambers. - Ejection Fraction is now calculated based on the simulated end-diastolic and end-systolic volumes of the left ventricle. - The example program has been updated to display a live, detailed summary of the heart's mechanical status, showing pressures, volumes, and valve states as they change over time.
This commit is contained in:
google-labs-jules[bot]
+170
-19
@@ -2,6 +2,16 @@
|
||||
#include <random>
|
||||
#include <algorithm>
|
||||
#include <sstream>
|
||||
#include <vector>
|
||||
#include <string>
|
||||
#include <cmath>
|
||||
#include <numeric>
|
||||
#include <deque>
|
||||
|
||||
// Helper for Gaussian function to model EKG waves
|
||||
static double gaussian(double x, double mu, double sigma) {
|
||||
return exp(-0.5 * pow((x - mu) / sigma, 2));
|
||||
}
|
||||
|
||||
// Helper function for random fluctuations
|
||||
static double getFluctuation(double stddev) {
|
||||
@@ -11,33 +21,174 @@ static double getFluctuation(double stddev) {
|
||||
return d(gen);
|
||||
}
|
||||
|
||||
Heart::Heart(int id) : Organ(id, "Heart"), heartRate(75.0), bloodPressureSystolic(120.0), bloodPressureDiastolic(80.0) {}
|
||||
Heart::Heart(int id, int numLeads)
|
||||
: Organ(id, "Heart"),
|
||||
heartRate(75.0),
|
||||
measuredHeartRate(75.0),
|
||||
numLeads(numLeads),
|
||||
ekgHistorySize(200),
|
||||
totalTime_s(0.0),
|
||||
cardiacCyclePosition_s(0.0),
|
||||
lastRPeakTime_s(-1.0),
|
||||
rPeakDetectedInCycle(false),
|
||||
ejectionFraction(0.55) {
|
||||
|
||||
// Initialize Chambers
|
||||
leftAtrium.name = "Left Atrium";
|
||||
rightAtrium.name = "Right Atrium";
|
||||
leftVentricle.name = "Left Ventricle";
|
||||
rightVentricle.name = "Right Ventricle";
|
||||
|
||||
// Initialize Valves
|
||||
mitralValve.name = "Mitral Valve";
|
||||
tricuspidValve.name = "Tricuspid Valve";
|
||||
aorticValve.name = "Aortic Valve";
|
||||
pulmonaryValve.name = "Pulmonary Valve";
|
||||
|
||||
const std::vector<std::string> allLeadNames = {
|
||||
"I", "II", "III", "aVR", "aVL", "aVF",
|
||||
"V1", "V2", "V3", "V4", "V5", "V6"
|
||||
};
|
||||
|
||||
int leadsToCreate = std::min((int)allLeadNames.size(), numLeads);
|
||||
for (int i = 0; i < leadsToCreate; ++i) {
|
||||
leadNames.push_back(allLeadNames[i]);
|
||||
ekgData[allLeadNames[i]] = std::deque<double>();
|
||||
}
|
||||
}
|
||||
|
||||
void Heart::update(double deltaTime_s) {
|
||||
const double baseline_hr = 75.0;
|
||||
const double baseline_bp_systolic = 120.0;
|
||||
const double baseline_bp_diastolic = 80.0;
|
||||
const double theta = 0.1; // Mean reversion speed
|
||||
const double hr_stddev = 0.1;
|
||||
const double bp_stddev = 0.1;
|
||||
// --- Electrical Simulation Update ---
|
||||
totalTime_s += deltaTime_s;
|
||||
heartRate += getFluctuation(0.01); // Slow variation in underlying rate
|
||||
heartRate = std::max(60.0, std::min(heartRate, 100.0));
|
||||
double cycleDuration_s = 60.0 / heartRate;
|
||||
|
||||
heartRate += theta * (baseline_hr - heartRate) * deltaTime_s + getFluctuation(hr_stddev * deltaTime_s);
|
||||
bloodPressureSystolic += theta * (baseline_bp_systolic - bloodPressureSystolic) * deltaTime_s + getFluctuation(bp_stddev * deltaTime_s);
|
||||
bloodPressureDiastolic += theta * (baseline_bp_diastolic - bloodPressureDiastolic) * deltaTime_s + getFluctuation(bp_stddev * deltaTime_s);
|
||||
double oldCyclePosition = cardiacCyclePosition_s;
|
||||
cardiacCyclePosition_s += deltaTime_s;
|
||||
|
||||
heartRate = std::clamp(heartRate, 60.0, 100.0);
|
||||
bloodPressureSystolic = std::clamp(bloodPressureSystolic, 90.0, 120.0);
|
||||
bloodPressureDiastolic = std::clamp(bloodPressureDiastolic, 60.0, 80.0);
|
||||
// R-peak detection for measured heart rate
|
||||
double rPeakCycleTime = 0.22 * cycleDuration_s;
|
||||
if (oldCyclePosition < rPeakCycleTime && cardiacCyclePosition_s >= rPeakCycleTime && !rPeakDetectedInCycle) {
|
||||
if (lastRPeakTime_s > 0) {
|
||||
measuredHeartRate = 60.0 / (totalTime_s - lastRPeakTime_s);
|
||||
}
|
||||
lastRPeakTime_s = totalTime_s;
|
||||
rPeakDetectedInCycle = true;
|
||||
}
|
||||
|
||||
if (cardiacCyclePosition_s > cycleDuration_s) {
|
||||
cardiacCyclePosition_s -= cycleDuration_s;
|
||||
rPeakDetectedInCycle = false;
|
||||
}
|
||||
|
||||
double timeInCycle = cardiacCyclePosition_s / cycleDuration_s;
|
||||
double baseVoltage = simulateEkgWaveform(timeInCycle);
|
||||
for(const auto& leadName : leadNames) {
|
||||
double leadModifier = 1.0 - (std::distance(leadNames.begin(), std::find(leadNames.begin(), leadNames.end(), leadName)) * 0.1);
|
||||
ekgData[leadName].push_front(baseVoltage * leadModifier);
|
||||
if (ekgData[leadName].size() > ekgHistorySize) ekgData[leadName].pop_back();
|
||||
}
|
||||
|
||||
// --- Mechanical Simulation Update ---
|
||||
// Define pressures outside the heart
|
||||
double venousPressure = 5.0; // CVP
|
||||
double pulmonaryArteryPressure = 20.0;
|
||||
double aorticPressure = getAorticPressure();
|
||||
|
||||
// Determine chamber states based on EKG-timed cycle
|
||||
if (timeInCycle >= 0.0 && timeInCycle < 0.15) { // Atrial Systole
|
||||
leftAtrium.state = ChamberState::SYSTOLE;
|
||||
rightAtrium.state = ChamberState::SYSTOLE;
|
||||
// End of ventricular filling, capture EDV
|
||||
if(oldCyclePosition < 0.15 && timeInCycle >= 0.15) {
|
||||
leftVentricle.endDiastolicVolume_mL = leftVentricle.volume_mL;
|
||||
}
|
||||
} else {
|
||||
leftAtrium.state = ChamberState::DIASTOLE;
|
||||
rightAtrium.state = ChamberState::DIASTOLE;
|
||||
}
|
||||
|
||||
if (timeInCycle >= 0.20 && timeInCycle < 0.5) { // Ventricular Systole
|
||||
leftVentricle.state = ChamberState::SYSTOLE;
|
||||
rightVentricle.state = ChamberState::SYSTOLE;
|
||||
// End of ventricular ejection, capture ESV and calculate EF
|
||||
if(oldCyclePosition < 0.5 && timeInCycle >= 0.5) {
|
||||
leftVentricle.endSystolicVolume_mL = leftVentricle.volume_mL;
|
||||
if (leftVentricle.endDiastolicVolume_mL > 0) {
|
||||
ejectionFraction = (leftVentricle.endDiastolicVolume_mL - leftVentricle.endSystolicVolume_mL) / leftVentricle.endDiastolicVolume_mL;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
leftVentricle.state = ChamberState::DIASTOLE;
|
||||
rightVentricle.state = ChamberState::DIASTOLE;
|
||||
}
|
||||
|
||||
// Update chamber pressures based on state (very simplified model)
|
||||
leftAtrium.pressure_mmHg = (leftAtrium.state == ChamberState::SYSTOLE) ? 10.0 : 5.0;
|
||||
rightAtrium.pressure_mmHg = (rightAtrium.state == ChamberState::SYSTOLE) ? 7.0 : 2.0;
|
||||
leftVentricle.pressure_mmHg = (leftVentricle.state == ChamberState::SYSTOLE) ? 125.0 * sin((timeInCycle - 0.2) / 0.3 * M_PI) : 5.0;
|
||||
rightVentricle.pressure_mmHg = (rightVentricle.state == ChamberState::SYSTOLE) ? 25.0 * sin((timeInCycle - 0.2) / 0.3 * M_PI) : 2.0;
|
||||
|
||||
// Update Valve Status
|
||||
tricuspidValve.status = (rightAtrium.pressure_mmHg > rightVentricle.pressure_mmHg) ? ValveStatus::OPEN : ValveStatus::CLOSED;
|
||||
mitralValve.status = (leftAtrium.pressure_mmHg > leftVentricle.pressure_mmHg) ? ValveStatus::OPEN : ValveStatus::CLOSED;
|
||||
pulmonaryValve.status = (rightVentricle.pressure_mmHg > pulmonaryArteryPressure) ? ValveStatus::OPEN : ValveStatus::CLOSED;
|
||||
aorticValve.status = (leftVentricle.pressure_mmHg > aorticPressure) ? ValveStatus::OPEN : ValveStatus::CLOSED;
|
||||
|
||||
// Update Chamber Volumes (simplified blood flow)
|
||||
double flowRate = 500.0 * deltaTime_s; // mL/s
|
||||
if (mitralValve.status == ValveStatus::OPEN) leftVentricle.volume_mL += flowRate;
|
||||
if (tricuspidValve.status == ValveStatus::OPEN) rightVentricle.volume_mL += flowRate;
|
||||
if (aorticValve.status == ValveStatus::OPEN) leftVentricle.volume_mL -= flowRate * 1.5;
|
||||
if (pulmonaryValve.status == ValveStatus::OPEN) rightVentricle.volume_mL -= flowRate * 1.5;
|
||||
|
||||
// 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));
|
||||
}
|
||||
|
||||
double Heart::simulateEkgWaveform(double timeInCycle) {
|
||||
double p_time = 0.1, q_time = 0.2, r_time = 0.22, s_time = 0.24, t_time = 0.4;
|
||||
double p_amp = 0.15, q_amp = -0.1, r_amp = 1.0, s_amp = -0.25, t_amp = 0.3;
|
||||
double p_sigma = 0.04, qrs_sigma = 0.02, t_sigma = 0.06;
|
||||
double voltage = 0.0;
|
||||
voltage += p_amp * gaussian(timeInCycle, p_time, p_sigma);
|
||||
voltage += q_amp * gaussian(timeInCycle, q_time, qrs_sigma);
|
||||
voltage += r_amp * gaussian(timeInCycle, r_time, qrs_sigma);
|
||||
voltage += s_amp * gaussian(timeInCycle, s_time, qrs_sigma);
|
||||
voltage += t_amp * gaussian(timeInCycle, t_time, t_sigma);
|
||||
return voltage;
|
||||
}
|
||||
|
||||
std::string Heart::getSummary() const {
|
||||
std::stringstream ss;
|
||||
ss << "Type: " << organType << " (ID: " << organId << ")\n"
|
||||
<< " Heart Rate: " << heartRate << " bpm\n"
|
||||
<< " Blood Pressure: " << bloodPressureSystolic << "/" << bloodPressureDiastolic << " mmHg";
|
||||
ss.precision(2);
|
||||
ss << std::fixed;
|
||||
ss << "--- Heart Summary ---\n"
|
||||
<< "Heart Rate (Measured): " << getHeartRate() << " bpm\n"
|
||||
<< "Ejection Fraction: " << getEjectionFraction() * 100.0 << "%\n"
|
||||
<< "Aortic Pressure: " << getAorticPressure() << " mmHg\n\n"
|
||||
<< "--- Chambers ---\n"
|
||||
<< " LV Volume: " << leftVentricle.volume_mL << " mL\n"
|
||||
<< " LV Pressure: " << leftVentricle.pressure_mmHg << " mmHg\n"
|
||||
<< " RV Volume: " << rightVentricle.volume_mL << " mL\n"
|
||||
<< " RV Pressure: " << rightVentricle.pressure_mmHg << " mmHg\n\n"
|
||||
<< "--- Valves ---\n"
|
||||
<< " Aortic Valve: " << (aorticValve.status == ValveStatus::OPEN ? "OPEN" : "CLOSED") << "\n"
|
||||
<< " Mitral Valve: " << (mitralValve.status == ValveStatus::OPEN ? "OPEN" : "CLOSED") << "\n";
|
||||
return ss.str();
|
||||
}
|
||||
|
||||
double Heart::getHeartRate() const { return heartRate; }
|
||||
double Heart::getBloodPressureSystolic() const { return bloodPressureSystolic; }
|
||||
double Heart::getBloodPressureDiastolic() const { return bloodPressureDiastolic; }
|
||||
double Heart::getHeartRate() const { return measuredHeartRate; }
|
||||
const std::map<std::string, std::deque<double>>& Heart::getEkgData() const { return ekgData; }
|
||||
double Heart::getEjectionFraction() const { return ejectionFraction; }
|
||||
double Heart::getAorticPressure() const {
|
||||
// Aortic pressure decays over time, but gets "pumped up" by ventricular ejection
|
||||
// This is a placeholder for a more complex arterial model.
|
||||
if(aorticValve.status == ValveStatus::OPEN) {
|
||||
return leftVentricle.pressure_mmHg;
|
||||
}
|
||||
// Simplified diastolic pressure decay
|
||||
return 80.0 + 40.0 * exp(-cardiacCyclePosition_s);
|
||||
}
|
||||
|
||||
+14
-3
@@ -15,16 +15,17 @@
|
||||
#include <memory>
|
||||
|
||||
/**
|
||||
* @brief Initializes a new patient with baseline vital signs and a standard set of organs.
|
||||
* @brief Initializes a new patient with a specific number of heart leads.
|
||||
* @param patientId The ID for the new patient.
|
||||
* @param numHeartLeads The number of EKG leads for the heart.
|
||||
* @return A Patient struct with default healthy values.
|
||||
*/
|
||||
Patient initializePatient(int patientId) {
|
||||
Patient initializePatient(int patientId, int numHeartLeads) {
|
||||
Patient patient;
|
||||
patient.patientId = patientId;
|
||||
|
||||
// Initialize default organs
|
||||
patient.organs.push_back(std::make_unique<Heart>(1));
|
||||
patient.organs.push_back(std::make_unique<Heart>(1, numHeartLeads));
|
||||
patient.organs.push_back(std::make_unique<Lungs>(2));
|
||||
patient.organs.push_back(std::make_unique<Brain>(3));
|
||||
patient.organs.push_back(std::make_unique<Liver>(4));
|
||||
@@ -41,6 +42,16 @@ Patient initializePatient(int patientId) {
|
||||
return patient;
|
||||
}
|
||||
|
||||
/**
|
||||
* @brief Initializes a new patient with baseline vital signs and a 12-lead heart.
|
||||
* @param patientId The ID for the new patient.
|
||||
* @return A Patient struct with default healthy values.
|
||||
*/
|
||||
Patient initializePatient(int patientId) {
|
||||
return initializePatient(patientId, 12);
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* @brief Updates the patient's state by updating the state of all their organs.
|
||||
* @param patient The patient to update.
|
||||
|
||||
Reference in New Issue
Block a user