Commit aac18f56 authored by Aaron Bray's avatar Aaron Bray
Browse files

Merge branch 'feature/baroreceptors' into '2.x'

Add Secondary Baroreceptor Effects

See merge request physiology/engine!25
parents ec0e0ca4 a7ce654a
......@@ -36,6 +36,7 @@
#include <fstream>
#include <stdexcept>
#include <algorithm>
#include <cmath>
#ifdef _MSC_VER
#include <direct.h>
......
......@@ -21,7 +21,6 @@
#include "properties/SEScalarMassPerVolume.h"
#include "properties/SEScalarAmountPerVolume.h"
#include <cmath>
#include <iterator>
DataTrack::DataTrack() : DataTrack(nullptr)
......
......@@ -14,7 +14,7 @@ protected:
/** @brief - Default system values to their homeostatic values */
virtual void Initialize() { SetUp(); }
virtual void SetUp() = 0;
virtual void SetUp() = 0; // Called after Initialize if stablizing, or after serialization if reading in a state
/// Notify systems that steady state has been achieved
virtual void AtSteadyState() {};
......
......@@ -322,6 +322,9 @@ void PBPulsePhysiology::Serialize(const NervousData& src, Nervous& dst)
dst.m_FeedbackActive = true;
dst.m_ArterialOxygenBaseline_mmHg = src.arterialoxygenbaseline_mmhg();
dst.m_ArterialCarbonDioxideBaseline_mmHg = src.arterialcarbondioxidebaseline_mmhg();
dst.m_BaroreceptorFeedbackStatus = src.baroreceptorfeedbackstatus();
dst.m_BaroreceptorActiveTime_s = src.baroreceptoractivetime_s();
dst.m_BaroreceptorMeanArterialPressureBaseline_mmHg = src.baroreceptormeanarterialpressurebaseline_mmhg();
}
NervousData* PBPulsePhysiology::Unload(const Nervous& src)
{
......@@ -334,6 +337,9 @@ void PBPulsePhysiology::Serialize(const Nervous& src, NervousData& dst)
PBPhysiology::Serialize(src, *dst.mutable_common());
dst.set_arterialoxygenbaseline_mmhg(src.m_ArterialOxygenBaseline_mmHg);
dst.set_arterialcarbondioxidebaseline_mmhg(src.m_ArterialCarbonDioxideBaseline_mmHg);
dst.set_baroreceptorfeedbackstatus(src.m_BaroreceptorFeedbackStatus);
dst.set_baroreceptoractivetime_s(src.m_BaroreceptorActiveTime_s);
dst.set_baroreceptormeanarterialpressurebaseline_mmhg(src.m_BaroreceptorMeanArterialPressureBaseline_mmHg);
}
void PBPulsePhysiology::Load(const RenalData& src, Renal& dst)
......
......@@ -2133,7 +2133,7 @@ void Cardiovascular::CalculatePleuralCavityVenousEffects()
{
//The left and right pleural pressures are likely to have large differences only due to a pneumothorax
//Pressure difference causes a mediastinum shift
double pleuralCavityPressureDiff_cmH2O = abs(m_leftPleuralCavity->GetPressure(PressureUnit::cmH2O) - m_rightPleuralCavity->GetPressure(PressureUnit::cmH2O));
double pleuralCavityPressureDiff_cmH2O = std::abs(m_leftPleuralCavity->GetPressure(PressureUnit::cmH2O) - m_rightPleuralCavity->GetPressure(PressureUnit::cmH2O));
double maxPressureDiff_cmH2O = 40.0;
double maxResistanceMultiplier = 10.0;
......
......@@ -71,6 +71,9 @@ void Nervous::Initialize()
// Set when feedback is turned on
m_ArterialOxygenBaseline_mmHg = 0;
m_ArterialCarbonDioxideBaseline_mmHg = 0;
m_BaroreceptorActiveTime_s = 0.0;
m_BaroreceptorFeedbackStatus = false;
m_BaroreceptorMeanArterialPressureBaseline_mmHg = 0;
}
//--------------------------------------------------------------------------------------------------
......@@ -83,15 +86,15 @@ void Nervous::Initialize()
void Nervous::SetUp()
{
m_dt_s = m_data.GetTimeStep().GetValue(TimeUnit::s);
m_normalizedGammaHeartRate = m_data.GetConfiguration().GetNormalizedHeartRateIntercept();
m_normalizedGammaElastance = m_data.GetConfiguration().GetNormalizedHeartElastanceIntercept();
m_normalizedGammaCompliance = m_data.GetConfiguration().GetNormalizedComplianceIntercept();
m_normalizedGammaResistance = m_data.GetConfiguration().GetNormalizedResistanceIntercept();
m_normalizedAlphaHeartRate = m_data.GetConfiguration().GetNormalizedHeartRateSympatheticSlope();
m_normalizedAlphaElastance = m_data.GetConfiguration().GetNormalizedHeartElastanceSympatheticSlope();
m_normalizedAlphaCompliance = m_data.GetConfiguration().GetNormalizedComplianceParasympatheticSlope();
m_normalizedAlphaResistance = m_data.GetConfiguration().GetNormalizedResistanceSympatheticSlope();
m_normalizedBetaHeartRate = m_data.GetConfiguration().GetNormalizedHeartRateParasympatheticSlope();
m_NormalizedGammaHeartRate = m_data.GetConfiguration().GetNormalizedHeartRateIntercept();
m_NormalizedGammaElastance = m_data.GetConfiguration().GetNormalizedHeartElastanceIntercept();
m_NormalizedGammaCompliance = m_data.GetConfiguration().GetNormalizedComplianceIntercept();
m_NormalizedGammaResistance = m_data.GetConfiguration().GetNormalizedResistanceIntercept();
m_NormalizedAlphaHeartRate = m_data.GetConfiguration().GetNormalizedHeartRateSympatheticSlope();
m_NormalizedAlphaElastance = m_data.GetConfiguration().GetNormalizedHeartElastanceSympatheticSlope();
m_NormalizedAlphaCompliance = m_data.GetConfiguration().GetNormalizedComplianceParasympatheticSlope();
m_NormalizedAlphaResistance = m_data.GetConfiguration().GetNormalizedResistanceSympatheticSlope();
m_NormalizedBetaHeartRate = m_data.GetConfiguration().GetNormalizedHeartRateParasympatheticSlope();
}
void Nervous::AtSteadyState()
......@@ -102,6 +105,7 @@ void Nervous::AtSteadyState()
// The set-points (Baselines) get reset at the end of each stabilization period.
m_ArterialOxygenBaseline_mmHg = m_data.GetBloodChemistry().GetArterialOxygenPressure(PressureUnit::mmHg);
m_ArterialCarbonDioxideBaseline_mmHg = m_data.GetBloodChemistry().GetArterialCarbonDioxidePressure(PressureUnit::mmHg);
m_BaroreceptorMeanArterialPressureBaseline_mmHg = m_data.GetPatient().GetMeanArterialPressureBaseline().GetValue(PressureUnit::mmHg);
// The baroreceptor scales need to be reset any time the baselines are reset.
GetBaroreceptorHeartRateScale().SetValue(1.0);
GetBaroreceptorHeartElastanceScale().SetValue(1.0);
......@@ -119,8 +123,11 @@ void Nervous::AtSteadyState()
//--------------------------------------------------------------------------------------------------
void Nervous::PreProcess()
{
BaroreceptorFeedback();
ChemoreceptorFeedback();
if (m_FeedbackActive)
{
BaroreceptorFeedback();
ChemoreceptorFeedback();
}
}
//--------------------------------------------------------------------------------------------------
......@@ -158,47 +165,80 @@ void Nervous::PostProcess()
/// These fractions are used to update the scaling parameters of heart rate, heart elastance, resistance and compliance
/// for each time step.
//--------------------------------------------------------------------------------------------------
/// \todo Add decompensation. Perhaps a reduction in the effect that is a function of blood volume below a threshold... and maybe time.
void Nervous::BaroreceptorFeedback()
{
if (!m_FeedbackActive)
return;
//First calculate the sympathetic and parasympathetic firing rates:
double nu = m_data.GetConfiguration().GetResponseSlope();
double meanArterialPressure_mmHg = m_data.GetCardiovascular().GetMeanArterialPressure(PressureUnit::mmHg);
//Adjusting the mean arterial pressure set-point to account for cardiovascular drug effects
double meanArterialPressureBaseline_mmHg = m_data.GetPatient().GetMeanArterialPressureBaseline(PressureUnit::mmHg) //m_MeanArterialPressureNoFeedbackBaseline_mmHg
double meanArterialPressureCombinedBaseline_mmHg = m_BaroreceptorMeanArterialPressureBaseline_mmHg
+ m_data.GetDrugs().GetMeanBloodPressureChange(PressureUnit::mmHg)
+ m_data.GetEnergy().GetExerciseMeanArterialPressureDelta(PressureUnit::mmHg);
double sympatheticFraction = 1.0 / (1.0 + pow(meanArterialPressure_mmHg / meanArterialPressureBaseline_mmHg, nu));
double parasympatheticFraction = 1.0 / (1.0 + pow(meanArterialPressure_mmHg / meanArterialPressureBaseline_mmHg, -nu));
+ m_data.GetEnergy().GetExerciseMeanArterialPressureDelta(PressureUnit::mmHg);
//Determine if barorectors are active and for how long
double pressureDeviation = meanArterialPressure_mmHg - meanArterialPressureCombinedBaseline_mmHg;
double pressureDeviationFraction = std::abs(pressureDeviation) / meanArterialPressureCombinedBaseline_mmHg;
if (pressureDeviationFraction > 0.05)
{
if (!m_BaroreceptorFeedbackStatus)
{
m_BaroreceptorActiveTime_s = 0.0;
m_BaroreceptorFeedbackStatus = true;
}
else
{
m_BaroreceptorActiveTime_s += m_dt_s;
//Need to change the threshold every 7 minutes based on over activation
if (m_BaroreceptorActiveTime_s > 420.0/*s*/)
{
//Reset the time
m_BaroreceptorActiveTime_s = 0.0;
//Adjust the threshold by up to 30%
m_BaroreceptorMeanArterialPressureBaseline_mmHg += 0.35 * pressureDeviation;
m_ss << "Baroreceptor MAP Baseline updated to " << m_BaroreceptorMeanArterialPressureBaseline_mmHg << " mmHg";
Info(m_ss);
}
}
}
else
{
m_BaroreceptorFeedbackStatus = false;
}
// \todo Dampen Baroreceptor response due to sedation
//double sedationDampeningEffect = m_data.GetDrugs().GetSedationLevel().GetValue() * 0.35; //maximum slope reduction of 35% - making that up right now
//update nu - the slope response of the firing rate
//nu *= (1 - sedationDampeningEffect); TODO: Add this back when I have data for validation
double sympatheticFraction = 1.0 / (1.0 + pow(meanArterialPressure_mmHg / meanArterialPressureCombinedBaseline_mmHg, nu));
double parasympatheticFraction = 1.0 / (1.0 + pow(meanArterialPressure_mmHg / meanArterialPressureCombinedBaseline_mmHg, -nu));
//Calculate the normalized change in heart rate
double normalizedHeartRate = GetBaroreceptorHeartRateScale().GetValue();
double tauHeartRate_s = m_data.GetConfiguration().GetHeartRateDistributedTimeDelay(TimeUnit::s);
double deltaNormalizedHeartRate = (1.0 / tauHeartRate_s)*(-normalizedHeartRate + m_normalizedAlphaHeartRate*sympatheticFraction - m_normalizedBetaHeartRate*parasympatheticFraction + m_normalizedGammaHeartRate)*m_dt_s;
double deltaNormalizedHeartRate = (1.0 / tauHeartRate_s)*(-normalizedHeartRate + m_NormalizedAlphaHeartRate*sympatheticFraction - m_NormalizedBetaHeartRate*parasympatheticFraction + m_NormalizedGammaHeartRate)*m_dt_s;
normalizedHeartRate += deltaNormalizedHeartRate;
GetBaroreceptorHeartRateScale().SetValue(normalizedHeartRate);
//Calculate the normalized change in heart elastance
double normalizedHeartElastance = GetBaroreceptorHeartElastanceScale().GetValue();
double tauElastance_s = m_data.GetConfiguration().GetHeartElastanceDistributedTimeDelay(TimeUnit::s);
double deltaNormalizedHeartElastance = (1.0 / tauElastance_s)*(-normalizedHeartElastance + m_normalizedAlphaElastance*sympatheticFraction + m_normalizedGammaElastance)*m_dt_s;
double deltaNormalizedHeartElastance = (1.0 / tauElastance_s)*(-normalizedHeartElastance + m_NormalizedAlphaElastance*sympatheticFraction + m_NormalizedGammaElastance)*m_dt_s;
normalizedHeartElastance += deltaNormalizedHeartElastance;
GetBaroreceptorHeartElastanceScale().SetValue(normalizedHeartElastance);
//Calculate the normalized change in flow resistance for any cardiovascular resistor
double normalizedResistance = GetBaroreceptorResistanceScale().GetValue();
double tauResistance_s = m_data.GetConfiguration().GetSystemicResistanceDistributedTimeDelay(TimeUnit::s);
double deltaNormalizedResistance = (1.0 / tauResistance_s)*(-normalizedResistance + m_normalizedAlphaResistance*sympatheticFraction + m_normalizedGammaResistance)*m_dt_s;
double deltaNormalizedResistance = (1.0 / tauResistance_s)*(-normalizedResistance + m_NormalizedAlphaResistance*sympatheticFraction + m_NormalizedGammaResistance)*m_dt_s;
normalizedResistance += deltaNormalizedResistance;
GetBaroreceptorResistanceScale().SetValue(normalizedResistance);
//Calculate the normalized change in flow compliance for any cardiovascular compliance
double normalizedCompliance = GetBaroreceptorComplianceScale().GetValue();
double tauCompliance_s = m_data.GetConfiguration().GetVenousComplianceDistributedTimeDelay(TimeUnit::s);
double deltaNormalizedCompliance = (1.0 / tauCompliance_s)*(-normalizedCompliance + m_normalizedAlphaCompliance*parasympatheticFraction + m_normalizedGammaCompliance)*m_dt_s;
double deltaNormalizedCompliance = (1.0 / tauCompliance_s)*(-normalizedCompliance + m_NormalizedAlphaCompliance*parasympatheticFraction + m_NormalizedGammaCompliance)*m_dt_s;
normalizedCompliance += deltaNormalizedCompliance;
GetBaroreceptorComplianceScale().SetValue(normalizedCompliance);
#ifdef VERBOSE
......@@ -258,9 +298,6 @@ void Nervous::CheckBrainStatus()
//--------------------------------------------------------------------------------------------------
void Nervous::ChemoreceptorFeedback()
{
if (!m_FeedbackActive)
return;
double normalized_pO2 = m_data.GetBloodChemistry().GetArterialOxygenPressure(PressureUnit::mmHg) / m_ArterialOxygenBaseline_mmHg;
double normalized_pCO2 = m_data.GetBloodChemistry().GetArterialCarbonDioxidePressure(PressureUnit::mmHg) / m_ArterialCarbonDioxideBaseline_mmHg;
......
......@@ -28,7 +28,7 @@ protected:
// Set pointers and other member variables common to both homeostatic initialization and loading a state
void SetUp();
void AtSteadyState();
void AtSteadyState();
void PreProcess();
void Process();
void PostProcess();
......@@ -38,21 +38,26 @@ protected:
void ChemoreceptorFeedback();
void SetPupilEffects();
// Serializable member variables (Set in Initialize and in schema)
// Serializable member variables (Set in Initialize and in schema)
double m_ArterialOxygenBaseline_mmHg;
double m_ArterialCarbonDioxideBaseline_mmHg;
//Baroreceptor threshold parameters
bool m_BaroreceptorFeedbackStatus;
double m_BaroreceptorActiveTime_s;
double m_BaroreceptorMeanArterialPressureBaseline_mmHg;
// Stateless member variable (Set in SetUp())
bool m_FeedbackActive;
bool m_FeedbackActive; // An overall flag used to turn off baro/chemo receptors during stabilization
double m_dt_s;
// Configuration fractions representing the amount by which each quantity can change due to baroreceptors;
double m_normalizedGammaHeartRate;
double m_normalizedGammaElastance;
double m_normalizedGammaCompliance;
double m_normalizedGammaResistance;
double m_normalizedAlphaHeartRate;
double m_normalizedAlphaElastance;
double m_normalizedAlphaCompliance;
double m_normalizedAlphaResistance;
double m_normalizedBetaHeartRate;
double m_NormalizedGammaHeartRate;
double m_NormalizedGammaElastance;
double m_NormalizedGammaCompliance;
double m_NormalizedGammaResistance;
double m_NormalizedAlphaHeartRate;
double m_NormalizedAlphaElastance;
double m_NormalizedAlphaCompliance;
double m_NormalizedAlphaResistance;
double m_NormalizedBetaHeartRate;
};
......@@ -16,8 +16,6 @@
#include "properties/SEScalarInversePressure.h"
#include "unsupported/Eigen/NonLinearOptimization"
#include <cmath>
#include <algorithm>
//#define VERBOSE
......@@ -315,7 +313,7 @@ void SaturationCalculator::CalculateCarbonMonoxideSpeciesDistribution(SELiquidCo
}
// Check to make sure hemoglobin was conserved
newTotalHb_mM = HbUnbound_mM + HbO2_mM + HbO2CO2_mM + HbCO2_mM + targetBoundCO_mM;
if (abs(newTotalHb_mM - totalHb_mM) > tolerance)
if (std::abs(newTotalHb_mM - totalHb_mM) > tolerance)
Warning("Hemoglobin not conserved during carbon monoxide species distribution calculation.");
// We can now set and balance for CO.
......@@ -430,12 +428,12 @@ void SaturationCalculator::CalculateBloodGasDistribution(SELiquidCompartment& cm
double newHb_mM = m_subHbQ->GetMolarity(AmountPerVolumeUnit::mmol_Per_L);
double newTotalHb_mM = newHbO2_mM + newHbCO2_mM + newHbO2CO2_mM + newHb_mM + newHbCO_mM;
double diffTotal = newTotalHb_mM - oldTotalHb_mM;
if (abs(diffTotal) > 1.0e-8)
if (std::abs(diffTotal) > 1.0e-8)
{
std::stringstream debugSS;
debugSS << "CalculateCarbonMonoxideSpeciesDistribution failed to conserve hemoglobin. Difference = ";
debugSS << diffTotal;
if (abs(diffTotal) > 1.0e-8)
if (std::abs(diffTotal) > 1.0e-8)
Warning(debugSS.str());
}
}
......
......@@ -95,6 +95,9 @@ message NervousData
cdm.NervousSystemData Common = 1;
double ArterialOxygenBaseline_mmHg = 2;
double ArterialCarbonDioxideBaseline_mmHg = 3;
bool BaroreceptorFeedbackStatus = 4;
double BaroreceptorActiveTime_s = 5;
double BaroreceptorMeanArterialPressureBaseline_mmHg = 6;
}
message RenalData
......
......@@ -54,9 +54,11 @@ patient/HemorrhageClass2FemoralArtery.json = ScenarioTest
patient/HemorrhageClass2InternalMultiple.json = ScenarioTest
patient/HemorrhageClass2InternalSpleen.json = ScenarioTest
patient/HemorrhageClass2NoFluid.json = ScenarioTest
patient/HemorrhageClass2NoFluidLongTerm.json = ScenarioTest
patient/HemorrhageClass2Saline.json = ScenarioTest
patient/HemorrhageClass3PackedRBC.json = ScenarioTest
patient/HemorrhageClass4NoFluid.json = ScenarioTest
patient/HemorrhageClass4NoFluidLongTerm.json = ScenarioTest
patient/HemorrhageCM1.json = ScenarioTest
patient/HemorrhageCM2.json = ScenarioTest
patient/LobarPneumoniaModerateBothLungs.json = ScenarioTest
......
422a6a97c6730b8ca2dcd948a961085b094eb7ae4468ef8ee406f07dede74c6581780095db7dc042bfe6ae9f06b9f26b97baf14895d915e514dacdd709c39ce8
\ No newline at end of file
eeadc52cd16b947e2b5f4388d77d46fea94b47a0a4be7e62f7ba7b70483902fa0998657e775b3c06cf6c77141536f5b9ce063de898669cec9901cf6eb00199a8
\ No newline at end of file
664dcda0f90cac2a4de4fcc979fb88fe1478510cce287d55a0b737527bcc0df2bf83b532e80eadeac9ef7505fba2576b8e47262b5adfc6e3c98717cf860f0821
\ No newline at end of file
933cf97a67d56e7eda770a48ab18b783487cdddbed5e5d77012d00ac2a367d11b5d266d96ca53898ffd0cf8a708238783de73ca07b39e9dc1d2471c560774a5a
\ No newline at end of file
e79a56ce649ddc481685bd1bd87851bb11ed404398ca29bed268ffd57a836e35b9120a4fc1bdcad206f3971c785896a6a16c76c0f566ae94cd1a699f934ce270
\ No newline at end of file
d03878a6ed1f2095f2121e2a9ac1da4abc6dad215025a31c6048f14078bf8e32a012b2b94c4a9930e11dfbcc4392ad2b9437c89d3f5c4fbfe28b1365240cc6d2
\ No newline at end of file
02c6142a425696a640a69425ef8f964d6eeb9c5eafb7dcc9cd66f997d762dc83879826391599b35a9cc05b2ca95ddaacbae951e5608fde794806530746c74b72
\ No newline at end of file
43b45e47e224b122e4f070c56051b700e0aadc597cc5b53de6a17ca0c25daed1d591dc4c95446fcf079960a876c12ea0da8665a80a17ca3f673c962319dae752
\ No newline at end of file
0c5fe8728b81a9f11204d7756084f6742b4a81e911983d218be28e6a9838d6196a4e66d7354390ff6504511c1a7f885c2234bb54635576bed9636f8f804a1ced
\ No newline at end of file
f869faf9ca3035d03d1caafb8599acaa1aeee740bd1ac6c3de7f5d0b403dd3ee82dccc5b90fbf8043172b44c766130dd76a9ca4e5ea987f2fcdc63959e810256
\ No newline at end of file
4cfdaf3f6037b7fbc499bb898c59e1d065e8fca52cf95797118f63d56276c0f079f08a158ec58e7fe9722614b7cb16b55f6e29b5e2e71fca8cdc40acd8e0fb5e
\ No newline at end of file
6f1a3681c314e35462660151074262381b15056c63f2a06163876dd6e436a4b062b0284c364ceb7dff785660a52ee2e4eec9fcd2736b1a934953cf9c2db6a0bb
\ No newline at end of file
7957e3c1bfe59c9f296a201b02961e3f057f451be465c350a9667e744a0732cb07572683c55e3a3bd920d03003ff429c0eb515d00830aadb74d887217aeb12df
\ No newline at end of file
ee10c82981bcb917c78de14ea28947c579de894b387f02a10950378d189fd2ad74a8ded3a9315b44808deed104b21609ba080f7ccba0ada709a7b006f68a7431
\ No newline at end of file
37a3ed3d98b8f6e9c7579a7fb7f979813c91796c0038a5e28abc9f50155e87be46c159944def9edffa3ae357978c611f9e85cfd05fd89c81b228a704d734b165
\ No newline at end of file
dd5eabc4083915b324b9258514264b52d830b5400788ca78f7517b4f581e6ca72f6649854918ba53bee616fe27f14de4fced24b4679c3ac1059d8d20e3040fae
\ No newline at end of file
bfdf218552e815815aed6004b1874d5bf4fa67a3984e8eef685325419611472ee8df993e198faa68a98674fc1e1d30e7c4d9a078b0b79132f80d5b374fbb57d8
\ No newline at end of file
3288d47cd9f41ad849e9f8a4c1d86a3657454a0e22353c3b6c4ba6d8edbba432482135d68f93d5451baef922dc0cb2540b3a9d59bc329f63d117a3fccc03aeb0
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment