Respiratory.cpp

/* Distributed under the Apache License, Version 2.0.
   See accompanying NOTICE file for details.*/

#include "stdafx.h"
#include "physiology/Respiratory.h"
#include "controller/Circuits.h"
#include "controller/Compartments.h"
#include "controller/Substances.h"
#include "PulseConfiguration.h"
// Conditions 
#include "engine/SEConditionManager.h"
#include "patient/conditions/SEChronicObstructivePulmonaryDisease.h"
#include "patient/conditions/SELobarPneumonia.h"
#include "patient/conditions/SEImpairedAlveolarExchange.h"
#include "patient/conditions/SEAcuteRespiratoryDistressSyndrome.h"
// Actions
#include "engine/SEActionManager.h"
#include "engine/SEPatientActionCollection.h"
#include "patient/actions/SEAcuteRespiratoryDistressSyndromeExacerbation.h"
#include "patient/actions/SEAirwayObstruction.h"
#include "patient/actions/SEAsthmaAttack.h"
#include "patient/actions/SEBronchoconstriction.h"
#include "patient/actions/SEChestOcclusiveDressing.h"
#include "patient/actions/SEChronicObstructivePulmonaryDiseaseExacerbation.h"
#include "patient/actions/SEConsciousRespiration.h"
#include "patient/actions/SEBreathHold.h"
#include "patient/actions/SEForcedExhale.h"
#include "patient/actions/SEForcedInhale.h"
#include "patient/actions/SEUseInhaler.h"
#include "patient/actions/SEDyspnea.h"
#include "patient/actions/SEIntubation.h"
#include "patient/actions/SELobarPneumoniaExacerbation.h"
#include "patient/actions/SEMechanicalVentilation.h"
#include "patient/actions/SENeedleDecompression.h"
#include "patient/actions/SERespiratoryFatigue.h"
#include "patient/actions/SESupplementalOxygen.h"
#include "patient/actions/SETensionPneumothorax.h"
// Assessments
#include "patient/assessments/SEPulmonaryFunctionTest.h"
// Dependent Systems
#include "system/physiology/SEBloodChemistrySystem.h"
#include "system/physiology/SEDrugSystem.h"
#include "system/physiology/SEEnergySystem.h"
// CDM
#include "patient/SEPatient.h"
#include "engine/SEEventManager.h"
#include "substance/SESubstance.h"
#include "substance/SESubstanceAerosolization.h"
#include "substance/SESubstanceConcentration.h"
#include "substance/SESubstanceFraction.h"
#include "substance/SESubstanceTransport.h"
#include "circuit/fluid/SEFluidCircuitCalculator.h"
#include "circuit/fluid/SEFluidCircuit.h"
#include "circuit/fluid/SEFluidCircuit.h"
#include "circuit/fluid/SEFluidCircuitNode.h"
#include "circuit/fluid/SEFluidCircuitPath.h"
#include "compartment/fluid/SEGasCompartmentGraph.h"
#include "compartment/fluid/SELiquidCompartmentGraph.h"
// Properties
#include "properties/SEScalar0To1.h"
#include "properties/SEScalarArea.h"
#include "properties/SEScalarFrequency.h"
#include "properties/SEScalarInversePressure.h"
#include "properties/SEScalarInverseVolume.h"
#include "properties/SEScalarLength.h"
#include "properties/SEScalarMass.h"
#include "properties/SEScalarMassPerVolume.h"
#include "properties/SEScalarNegative1To1.h"
#include "properties/SEScalarPower.h"
#include "properties/SEScalarPressure.h"
#include "properties/SEScalarPressurePerVolume.h"
#include "properties/SEScalarPressureTimePerVolume.h"
#include "properties/SEScalarTime.h"
#include "properties/SEScalarVolume.h"
#include "properties/SEScalarVolumePerTime.h"
#include "properties/SEScalarVolumePerPressure.h"
#include "properties/SEFunctionVolumeVsTime.h"
#include "properties/SERunningAverage.h"
// Data Track
#include "utils/DataTrack.h"

#ifdef _MSC_VER 
#pragma warning( disable : 4305 4244 )  // Disable some warning messages
#endif

//Flag for data tracks
//Should be commented out, unless debugging
//#define DEBUG

//Flag for setting things constant to test
//Should be commented out, unless debugging/tuning
//#define TUNING

Respiratory::Respiratory(PulseController& data) : SERespiratorySystem(data.GetLogger()), m_data(data)
{
  m_BloodPHRunningAverage = new SERunningAverage();
  m_ArterialO2RunningAverage_mmHg = new SERunningAverage();
  m_ArterialCO2RunningAverage_mmHg = new SERunningAverage();

  m_Calculator = new SEFluidCircuitCalculator(VolumePerPressureUnit::L_Per_cmH2O, VolumePerTimeUnit::L_Per_s, PressureTimeSquaredPerVolumeUnit::cmH2O_s2_Per_L, PressureUnit::cmH2O, VolumeUnit::L, PressureTimePerVolumeUnit::cmH2O_s_Per_L, GetLogger());
  m_GasTransporter = new SEGasTransporter(VolumePerTimeUnit::L_Per_s, VolumeUnit::L, VolumeUnit::L, GetLogger());
  m_AerosolTransporter = new SELiquidTransporter(VolumePerTimeUnit::mL_Per_s, VolumeUnit::mL, MassUnit::ug, MassPerVolumeUnit::ug_Per_mL, GetLogger());
  Clear();
}

Respiratory::~Respiratory()
{
  Clear();
  delete m_BloodPHRunningAverage;
  delete m_ArterialO2RunningAverage_mmHg;
  delete m_ArterialCO2RunningAverage_mmHg;

  delete m_Calculator;
  delete m_GasTransporter;
  delete m_AerosolTransporter;
}

void Respiratory::Clear()
{
  SERespiratorySystem::Clear();
  m_PatientActions = nullptr;

  m_Environment = nullptr;
  m_AerosolMouth = nullptr;
  m_AerosolCarina = nullptr;
  m_AerosolLeftAnatomicDeadSpace = nullptr;
  m_AerosolLeftAlveolarDeadSpace = nullptr;
  m_AerosolLeftAlveoli = nullptr;
  m_AerosolRightAnatomicDeadSpace = nullptr;
  m_AerosolRightAlveolarDeadSpace = nullptr;
  m_AerosolRightAlveoli = nullptr;
  m_Lungs = nullptr;
  m_LeftLungExtravascular = nullptr;
  m_RightLungExtravascular = nullptr;
  m_Carina = nullptr;
  m_AortaO2 = nullptr;
  m_AortaCO2 = nullptr;
  m_MechanicalVentilatorConnection = nullptr;
  m_MechanicalVentilatorAerosolConnection = nullptr;
  m_PleuralCavity = nullptr;

  m_RespiratoryCircuit = nullptr;

  m_Mouth = nullptr;
  m_LeftAlveoli = nullptr;
  m_LeftAnatomicDeadSpace = nullptr;
  m_LeftAlveolarDeadSpace = nullptr;
  m_LeftPleural = nullptr;
  m_RespiratoryMuscle = nullptr;
  m_RightAlveoli = nullptr;
  m_RightAnatomicDeadSpace = nullptr;
  m_RightAlveolarDeadSpace = nullptr;
  m_RightPleural = nullptr;

  m_CarinaToLeftAnatomicDeadSpace = nullptr;
  m_CarinaToRightAnatomicDeadSpace = nullptr;
  m_LeftAnatomicDeadSpaceToLeftAlveolarDeadSpace = nullptr;
  m_RightAnatomicDeadSpaceToRightAlveolarDeadSpace = nullptr;
  m_LeftAlveolarDeadSpaceToLeftAlveoli = nullptr;
  m_RightAlveolarDeadSpaceToRightAlveoli = nullptr;
  m_RightPleuralToRespiratoryMuscle = nullptr;
  m_LeftPleuralToRespiratoryMuscle = nullptr;
  m_DriverPressurePath = nullptr;
  m_MouthToCarina = nullptr;
  m_MouthToStomach = nullptr;
  m_EnvironmentToLeftChestLeak = nullptr;
  m_EnvironmentToRightChestLeak = nullptr;
  m_LeftAlveoliLeakToLeftPleural = nullptr;
  m_RightAlveoliLeakToRightPleural = nullptr;
  m_LeftPleuralToEnvironment = nullptr;
  m_RightPleuralToEnvironment = nullptr;
  m_RightAlveoliToRightPleuralConnection = nullptr;
  m_LeftAlveoliToLeftPleuralConnection = nullptr;
  m_RightPulmonaryCapillary = nullptr;
  m_LeftPulmonaryCapillary = nullptr;
  m_ConnectionToMouth = nullptr;
  m_GroundToConnection = nullptr;

  m_BloodPHRunningAverage->Clear();
  m_ArterialO2RunningAverage_mmHg->Clear();
  m_ArterialCO2RunningAverage_mmHg->Clear();
}

//--------------------------------------------------------------------------------------------------
/// \brief
/// Initializes system properties to valid homeostatic values.
//--------------------------------------------------------------------------------------------------
void Respiratory::Initialize()
{
  PulseSystem::Initialize();

  //Vital signs
  m_NotBreathing = false;
  m_TopBreathTotalVolume_L = 0.0;
  m_TopBreathAlveoliVolume_L = 0.0;
  m_TopBreathPleuralVolume_L = 0.0;
  m_TopBreathPleuralPressure_cmH2O = 0.0;
  m_TopBreathAlveoliPressure_cmH2O = 0.0;
  m_TopBreathDriverPressure_cmH2O = 0.0;
  m_LastCardiacCycleBloodPH = 7.4;
  m_TopCarinaO2 = 0.2;
  m_TopBreathElapsedTime_min = 0.0;
  m_BottomBreathElapsedTime_min = 0.0;
  m_BottomBreathTotalVolume_L = 0.0;
  m_BottomBreathAlveoliVolume_L = m_BottomBreathTotalVolume_L;
  m_BottomBreathPleuralVolume_L = 0.0;
  m_BottomBreathPleuralPressure_cmH2O = 0.0;
  m_BottomBreathAlveoliPressure_cmH2O = 0.0;
  m_BottomBreathDriverPressure_cmH2O = 0.0;
  m_PeakAlveolarPressure_cmH2O = 0.0;
  m_MaximalAlveolarPressure_cmH2O = 0.0;

  //Driver
  m_MaxDriverPressure_cmH2O = -50.0;
  m_ElapsedBreathingCycleTime_min = 0.0;
  m_TopBreathElapsedTime_min = 0.0;
  m_BreathingCycle = false;
  m_BreathingCycleTime_s = 0.0;
  m_VentilationFrequency_Per_min = m_data.GetCurrentPatient().GetRespirationRateBaseline(FrequencyUnit::Per_min);
  m_DriverPressure_cmH2O = 0.0;
  m_DriverPressureMin_cmH2O = 0.0;
  m_VentilationToTidalVolumeSlope = 30.0;
  //The peak driver pressure is the pressure above the default pressure
  m_PeakRespiratoryDrivePressure_cmH2O = 0.0;
  m_ArterialO2PartialPressure_mmHg = m_AortaO2->GetPartialPressure(PressureUnit::mmHg);
  m_ArterialCO2PartialPressure_mmHg = m_AortaCO2->GetPartialPressure(PressureUnit::mmHg);
  m_PreviousTargetAlveolarVentilation_L_Per_min = m_data.GetCurrentPatient().GetTidalVolumeBaseline(VolumeUnit::L) * m_VentilationFrequency_Per_min;
  m_AverageLocalTissueBronchodilationEffects = 0.0;

  m_IERatioScaleFactor = 1.0;

  // Conscious Respiration
  m_ConsciousRespirationPeriod_s = 0.0;
  m_ConsciousRespirationRemainingPeriod_s = 0.0;
  m_ExpiratoryReserveVolumeFraction = -1.0;
  m_InspiratoryCapacityFraction = -1.0;
  m_ConsciousStartPressure_cmH2O = 0.0;
  m_ConsciousEndPressure_cmH2O = 0.0;
  m_ConsciousBreathing = false;

  //System data
  double TidalVolume_L = m_data.GetCurrentPatient().GetTidalVolumeBaseline(VolumeUnit::L);
  double RespirationRate_Per_min = m_data.GetCurrentPatient().GetRespirationRateBaseline(FrequencyUnit::Per_min);
  GetTidalVolume().SetValue(TidalVolume_L, VolumeUnit::L);
  GetRespirationRate().SetValue(RespirationRate_Per_min, FrequencyUnit::Per_min);
  GetInspiratoryExpiratoryRatio().SetValue(0.5);
  GetCarricoIndex().SetValue(452.0, PressureUnit::mmHg);

  double AnatomicDeadSpace_L = m_LeftAnatomicDeadSpace->GetVolumeBaseline(VolumeUnit::L) + m_RightAnatomicDeadSpace->GetVolumeBaseline(VolumeUnit::L);
  double AlveolarDeadSpace_L = 0.0;
  GetAnatomicDeadSpace().SetValue(AnatomicDeadSpace_L, VolumeUnit::L);
  GetAlveolarDeadSpace().SetValue(AlveolarDeadSpace_L, VolumeUnit::L);
  GetPhysiologicDeadSpace().SetValue(AnatomicDeadSpace_L + AlveolarDeadSpace_L, VolumeUnit::L);

  GetTotalAlveolarVentilation().SetValue(RespirationRate_Per_min * TidalVolume_L, VolumePerTimeUnit::L_Per_min);
  GetTotalPulmonaryVentilation().SetValue(RespirationRate_Per_min * TidalVolume_L, VolumePerTimeUnit::L_Per_min);
  GetTotalDeadSpaceVentilation().SetValue(AnatomicDeadSpace_L * RespirationRate_Per_min, VolumePerTimeUnit::L_Per_min);
  GetSpecificVentilation().SetValue(0.21);

  GetEndTidalCarbonDioxideFraction().SetValue(0.0827);
  GetEndTidalCarbonDioxidePressure().SetValue(0.0, PressureUnit::mmHg);
  GetEndTidalOxygenFraction().SetValue(0.0);
  GetEndTidalOxygenPressure().SetValue(0.0, PressureUnit::mmHg);

  GetIntrapleuralPressure().SetValue(0.0, PressureUnit::cmH2O);
  GetTransrespiratoryPressure().SetValue(0.0, PressureUnit::cmH2O);
  GetTransairwayPressure().SetValue(0.0, PressureUnit::cmH2O);
  GetTranspulmonaryPressure().SetValue(0.0, PressureUnit::cmH2O);
  GetTransalveolarPressure().SetValue(0.0, PressureUnit::cmH2O);
  GetTransthoracicPressure().SetValue(0.0, PressureUnit::cmH2O);
  GetTransChestWallPressure().SetValue(0.0, PressureUnit::cmH2O);

  GetPulmonaryCompliance().SetValue(0.1, VolumePerPressureUnit::L_Per_cmH2O);
  GetLungCompliance().SetValue(0.1, VolumePerPressureUnit::L_Per_cmH2O);
  GetChestWallCompliance().SetValue(0.2, VolumePerPressureUnit::L_Per_cmH2O);
  GetPulmonaryElastance().SetValue(1.0 / 0.1, PressurePerVolumeUnit::cmH2O_Per_L);

  // Muscle Pressure Waveform
  m_InspiratoryRiseFraction = 0;
  m_InspiratoryHoldFraction = 0;
  m_InspiratoryReleaseFraction = 0;
  m_InspiratoryToExpiratoryPauseFraction = 0;
  m_ExpiratoryRiseFraction = 0;
  m_ExpiratoryHoldFraction = 0;
  m_ExpiratoryReleaseFraction = 0;
  m_ResidueFraction = 0;

  //Get the fluid mechanics to a good starting point
  TuneCircuit();
}

//--------------------------------------------------------------------------------------------------
/// \brief
/// Initializes parameters for Respiratory Class
///
///  \details
///   Initializes member variables and system level values on the common data model.
//--------------------------------------------------------------------------------------------------
void Respiratory::SetUp()
{
  //Time Step
  m_dt_s = m_data.GetTimeStep().GetValue(TimeUnit::s);
  m_dt_min = m_data.GetTimeStep().GetValue(TimeUnit::min);
  //Patient
  m_PatientActions = &m_data.GetActions().GetPatientActions();
  //Configuration parameters
  m_DefaultOpenResistance_cmH2O_s_Per_L = m_data.GetConfiguration().GetDefaultOpenFlowResistance(PressureTimePerVolumeUnit::cmH2O_s_Per_L);
  m_DefaultClosedResistance_cmH2O_s_Per_L = m_data.GetConfiguration().GetDefaultClosedFlowResistance(PressureTimePerVolumeUnit::cmH2O_s_Per_L);
  m_RespOpenResistance_cmH2O_s_Per_L = m_data.GetConfiguration().GetRespiratoryOpenResistance(PressureTimePerVolumeUnit::cmH2O_s_Per_L);
  m_RespClosedResistance_cmH2O_s_Per_L = m_data.GetConfiguration().GetRespiratoryClosedResistance(PressureTimePerVolumeUnit::cmH2O_s_Per_L);
  m_PeripheralControlGainConstant = m_data.GetConfiguration().GetPeripheralVentilatoryControllerGain();
  m_CentralControlGainConstant = m_data.GetConfiguration().GetCentralVentilatoryControllerGain();
  m_VentilationTidalVolumeIntercept = m_data.GetConfiguration().GetVentilationTidalVolumeIntercept(VolumeUnit::L);
  m_VentilatoryOcclusionPressure_cmH2O = m_data.GetConfiguration().GetVentilatoryOcclusionPressure(PressureUnit::cmH2O); //This increases the absolute max driver pressure
  m_PleuralComplianceSensitivity_Per_L = m_data.GetConfiguration().GetPleuralComplianceSensitivity(InverseVolumeUnit::Inverse_L);
  m_MinimumAllowableTidalVolume_L = m_data.GetConfiguration().GetMinimumAllowableTidalVolume(VolumeUnit::L);
  m_MinimumAllowableInpiratoryAndExpiratoryPeriod_s = m_data.GetConfiguration().GetMinimumAllowableInpiratoryAndExpiratoryPeriod(TimeUnit::s);
  //Compartments
  m_Environment = m_data.GetCompartments().GetGasCompartment(pulse::EnvironmentCompartment::Ambient);
  m_AerosolMouth = m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::Mouth);
  m_AerosolCarina = m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::Carina);
  m_AerosolLeftAnatomicDeadSpace = m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::LeftAnatomicDeadSpace);
  m_AerosolLeftAlveolarDeadSpace = m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::LeftAlveolarDeadSpace);
  m_AerosolLeftAlveoli = m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::LeftAlveoli);
  m_AerosolRightAnatomicDeadSpace = m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::RightAnatomicDeadSpace);
  m_AerosolRightAlveolarDeadSpace = m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::RightAlveolarDeadSpace);
  m_AerosolRightAlveoli = m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::RightAlveoli);
  m_Lungs = m_data.GetCompartments().GetGasCompartment(pulse::PulmonaryCompartment::Lungs);
  m_LeftLung = m_data.GetCompartments().GetGasCompartment(pulse::PulmonaryCompartment::LeftLung);
  m_RightLung = m_data.GetCompartments().GetGasCompartment(pulse::PulmonaryCompartment::RightLung);
  m_PleuralCavity = m_data.GetCompartments().GetGasCompartment(pulse::PulmonaryCompartment::PleuralCavity);
  m_LeftLungExtravascular = m_data.GetCompartments().GetLiquidCompartment(pulse::ExtravascularCompartment::LeftLungIntracellular);
  m_RightLungExtravascular = m_data.GetCompartments().GetLiquidCompartment(pulse::ExtravascularCompartment::RightLungIntracellular);
  m_Carina = m_data.GetCompartments().GetGasCompartment(pulse::PulmonaryCompartment::Carina);
  m_CarinaO2 = m_Carina->GetSubstanceQuantity(m_data.GetSubstances().GetO2());
  SELiquidCompartment* Aorta = m_data.GetCompartments().GetLiquidCompartment(pulse::VascularCompartment::Aorta);
  m_AortaO2 = Aorta->GetSubstanceQuantity(m_data.GetSubstances().GetO2());
  m_AortaCO2 = Aorta->GetSubstanceQuantity(m_data.GetSubstances().GetCO2());
  m_LeftAlveoliO2 = m_data.GetCompartments().GetGasCompartment(pulse::PulmonaryCompartment::LeftAlveoli)->GetSubstanceQuantity(m_data.GetSubstances().GetO2());
  m_RightAlveoliO2 = m_data.GetCompartments().GetGasCompartment(pulse::PulmonaryCompartment::RightAlveoli)->GetSubstanceQuantity(m_data.GetSubstances().GetO2());
  m_MechanicalVentilatorConnection = m_data.GetCompartments().GetGasCompartment(pulse::MechanicalVentilatorCompartment::Connection);
  m_MechanicalVentilatorAerosolConnection = m_data.GetCompartments().GetLiquidCompartment(pulse::MechanicalVentilatorCompartment::Connection);
  // Compartments we will process aerosol effects on
  m_AerosolEffects.clear();
  m_AerosolEffects.push_back(m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::Carina));
  m_AerosolEffects.push_back(m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::LeftAlveoli));
  m_AerosolEffects.push_back(m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::LeftAnatomicDeadSpace));
  m_AerosolEffects.push_back(m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::LeftAlveolarDeadSpace));
  m_AerosolEffects.push_back(m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::RightAlveoli));
  m_AerosolEffects.push_back(m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::RightAnatomicDeadSpace));
  m_AerosolEffects.push_back(m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::RightAlveolarDeadSpace));
  //Circuits
  m_RespiratoryCircuit = &m_data.GetCircuits().GetRespiratoryCircuit();
  //Nodes
  m_Mouth = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::Mouth);
  m_LeftAlveoli = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::LeftAlveoli);
  m_LeftAnatomicDeadSpace = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::LeftAnatomicDeadSpace);
  m_LeftAlveolarDeadSpace = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::LeftAlveolarDeadSpace);
  m_LeftPleural = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::LeftPleural);
  m_RespiratoryMuscle = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::RespiratoryMuscle);
  m_RightAlveoli = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::RightAlveoli);
  m_RightAnatomicDeadSpace = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::RightAnatomicDeadSpace);
  m_RightAlveolarDeadSpace = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::RightAlveolarDeadSpace);
  m_RightPleural = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::RightPleural);
  m_Ambient = m_RespiratoryCircuit->GetNode(pulse::EnvironmentNode::Ambient);
  m_Stomach = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::Stomach);
  //Paths
  m_CarinaToLeftAnatomicDeadSpace = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::CarinaToLeftAnatomicDeadSpace);
  m_CarinaToRightAnatomicDeadSpace = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::CarinaToRightAnatomicDeadSpace);
  m_LeftAnatomicDeadSpaceToLeftAlveolarDeadSpace = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::LeftAnatomicDeadSpaceToLeftAlveolarDeadSpace);
  m_RightAnatomicDeadSpaceToRightAlveolarDeadSpace = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::RightAnatomicDeadSpaceToRightAlveolarDeadSpace);
  m_LeftAlveolarDeadSpaceToLeftAlveoli = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::LeftAlveolarDeadSpaceToLeftAlveoli);
  m_RightAlveolarDeadSpaceToRightAlveoli = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::RightAlveolarDeadSpaceToRightAlveoli);
  m_RightPleuralToRespiratoryMuscle = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::RightPleuralToRespiratoryMuscle);
  m_LeftPleuralToRespiratoryMuscle = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::LeftPleuralToRespiratoryMuscle);
  m_DriverPressurePath = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::EnvironmentToRespiratoryMuscle);
  m_MouthToCarina = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::MouthToCarina);
  m_MouthToStomach = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::MouthToStomach);
  m_EnvironmentToLeftChestLeak = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::EnvironmentToLeftChestLeak);
  m_EnvironmentToRightChestLeak = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::EnvironmentToRightChestLeak);
  m_LeftAlveoliLeakToLeftPleural = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::LeftAlveoliLeakToLeftPleural);
  m_RightAlveoliLeakToRightPleural = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::RightAlveoliLeakToRightPleural);
  m_LeftPleuralToEnvironment = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::LeftPleuralToEnvironment);
  m_RightPleuralToEnvironment = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::RightPleuralToEnvironment);
  m_RightAlveoliToRightPleuralConnection = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::RightAlveoliToRightPleuralConnection);
  m_LeftAlveoliToLeftPleuralConnection = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::LeftAlveoliToLeftPleuralConnection);
  m_ConnectionToMouth = m_data.GetCircuits().GetRespiratoryAndMechanicalVentilatorCircuit().GetPath(pulse::MechanicalVentilatorPath::ConnectionToMouth);
  m_GroundToConnection = m_data.GetCircuits().GetRespiratoryAndMechanicalVentilatorCircuit().GetPath(pulse::MechanicalVentilatorPath::GroundToConnection);

  /// \todo figure out how to modify these resistances without getting the cv circuit - maybe add a parameter, like baroreceptors does
  m_RightPulmonaryCapillary = m_data.GetCircuits().GetCardiovascularCircuit().GetPath(pulse::CardiovascularPath::RightPulmonaryCapillariesToRightPulmonaryVeins);
  m_LeftPulmonaryCapillary = m_data.GetCircuits().GetCardiovascularCircuit().GetPath(pulse::CardiovascularPath::LeftPulmonaryCapillariesToLeftPulmonaryVeins);
}

//--------------------------------------------------------------------------------------------------
/// \brief
/// Respiratory system at steady state
///
/// \details
/// Initializes respiratory conditions if any are present.
///  <UL>
///   <LI>COPD</LI>
///   <LI>Lobar Pneumonia</LI>
///   <LI>ImpairedAlveolarExchange</LI>
///  </UL>
///
//--------------------------------------------------------------------------------------------------
void Respiratory::AtSteadyState()
{
  /// \todo figure out how to save an initial healthy/baseline patient definition and another current patient definition
  //Experimentally determined  
  double tidalVolumeBaseline_L = GetTidalVolume(VolumeUnit::L);
  //Calculated
  double totalLungCapacity_L = m_data.GetCurrentPatient().GetTotalLungCapacity(VolumeUnit::L);
  double functionalResidualCapacity_L = m_data.GetCurrentPatient().GetFunctionalResidualCapacity(VolumeUnit::L);
  double inspiratoryReserveVolume_L = totalLungCapacity_L - functionalResidualCapacity_L - tidalVolumeBaseline_L;

  std::string typeString = "Initial Stabilization Homeostasis: ";
  if (m_data.GetState() == EngineState::AtSecondaryStableState)
    typeString = "Secondary Stabilization Homeostasis: ";

  std::stringstream ss;
  ss << typeString << "Patient tidal volume = " << tidalVolumeBaseline_L << " L.";
  Info(ss);
  ss << typeString << "Patient inspiratory reserve volume = " << inspiratoryReserveVolume_L << " L.";
  Info(ss);

  if (m_data.GetState() == EngineState::AtInitialStableState)
  {
    //At Resting State
    //Note: Respiratory conditions are applied each timestep to handle combined effects properly

    //Update healthy patient volumes
    //Tidal Volume is experimentatlly determined and is dependent on other patient settings
    //Inspiratory Reserve Volume is calcualted from Tidal Volume
    m_data.GetCurrentPatient().GetTidalVolumeBaseline().SetValue(tidalVolumeBaseline_L, VolumeUnit::L);
    m_data.GetCurrentPatient().GetInspiratoryReserveVolume().SetValue(inspiratoryReserveVolume_L, VolumeUnit::L);
  }
}

//--------------------------------------------------------------------------------------------------
/// \brief
/// Respiratory system preprocess function
///
/// \details
/// Calculates the muscle pressure source that drives the respiratory system.
/// Updates airway resistance to account for changes arising from factors 
/// like drugs and respiratory insults and interventions. 
//// Updates the chest wall variable compliance. Handles all respiratory
/// insults and actions.
//--------------------------------------------------------------------------------------------------
void Respiratory::PreProcess()
{
  CalculateWork();
  CalculateFatigue();

  UpdateChestWallCompliances();
  UpdateVolumes();
  UpdateResistances();
  UpdateAlveolarCompliances();
  UpdateInspiratoryExpiratoryRatio();
  UpdateDiffusion();
  UpdatePulmonaryCapillary();

  ProcessAerosolSubstances();
  Pneumothorax();  
  ConsciousRespiration();

  MechanicalVentilation();
  SupplementalOxygen();

  RespiratoryDriver();
}

//--------------------------------------------------------------------------------------------------
/// \brief
/// Respiratory system process function
///
/// \details
/// Ensures the selection of the respiratory system with or without the anesthesia machine. 
///  Handles the integration of the anesthesia machine to the respiratory system when the anesthesia machine is turned on.
/// The integration of the anesthesia machine to the respiratory system is handled at run time by constructing a combined circuit of 
/// the respiratory and anesthesia machine.  
/// Handles lung volume changes during alveolar gas transfer. 
/// Calculates physiological parameters such as respiration rate, tidal volume and others that belonging to the respiratory system.
//--------------------------------------------------------------------------------------------------
void Respiratory::Process()
{
  // Respiration circuit changes based on if Anesthesia Machine is on or off
  // When dynamic intercircuit connections work, we can stash off the respiration circuit in a member variable
  SEFluidCircuit& RespirationCircuit = m_data.GetCircuits().GetActiveRespiratoryCircuit();
  
  // Calc the circuits
  m_Calculator->Process(RespirationCircuit, m_dt_s);
  
  //ModifyPleuralVolume();
  SEGasCompartmentGraph& RespirationGraph = m_data.GetCompartments().GetActiveRespiratoryGraph();
  SELiquidCompartmentGraph& AerosolGraph = m_data.GetCompartments().GetActiveAerosolGraph();
  
  // Transport substances
  m_GasTransporter->Transport(RespirationGraph, m_dt_s);
  if(m_AerosolMouth->HasSubstanceQuantities())
    m_AerosolTransporter->Transport(AerosolGraph, m_dt_s);
  
  //Update system data
  CalculateVitalSigns();

#ifdef DEBUG
  Debugging(RespirationCircuit);
#endif  
}

//--------------------------------------------------------------------------------------------------
/// \brief
/// Respiratory system postprocess function
///
/// \details
/// Updates the current values of the gas volume fraction and gas volumes for the nodes in the respiratory circuit 
/// or the nodes in the combined (respiratory + anesthesia machine) circuit when the anesthesia machine is turned on.
//--------------------------------------------------------------------------------------------------
void Respiratory::PostProcess()
{
  // Respiration circuit changes based on if Anesthesia Machine is on or off
  // When dynamic intercircuit connections work, we can stash off the respiration circuit in a member variable
  SEFluidCircuit& RespirationCircuit = m_data.GetCircuits().GetActiveRespiratoryCircuit();  
  m_Calculator->PostProcess(RespirationCircuit);  
}

//--------------------------------------------------------------------------------------------------
/// \brief
/// Compute deposited mass, update localized PD effects 
///
/// \details
/// For each aerosol get the SIDE coefficient to determine deposited mass in each respiratory compartment. 
/// Adjust the resistances between compartments as a function of deposited mass to reach validated data.  
/// Liquid and solid aerosols are handeled here. 
//--------------------------------------------------------------------------------------------------
void Respiratory::ProcessAerosolSubstances()
{
  m_AverageLocalTissueBronchodilationEffects = 0.0; //No effect

  size_t numAerosols = m_AerosolMouth->GetSubstanceQuantities().size();
  if (numAerosols == 0)
    return;

  double inflammationCoefficient;

  // For this time step
  double mouthDepositied_ug = 0;
  double carinaDepositied_ug = 0;
  double leftDeadSpaceDepositied_ug = 0;
  double leftAlveoliDepositied_ug = 0;
  double rightDeadSpaceDepositied_ug = 0;
  double rightAlveoliDepositied_ug = 0;

  // Total amount deposited (including this time step)
  double mouthTotalDepositied_ug = 0;
  double carinaTotalDepositied_ug = 0;
  double leftDeadSpaceTotalDepositied_ug = 0;
  double leftAlveoliTotalDepositied_ug = 0;
  double rightDeadSpaceTotalDepositied_ug = 0;
  double rightAlveoliTotalDepositied_ug = 0;

  // Resistance Modifier Sum
  double mouthResistanceModifier=1;
  double carinaResistanceModifier=1;
  double leftDeadSpaceResistanceModifier=1;
  double leftAlveoliResistanceModifier=1;
  double rightDeadSpaceResistanceModifier=1;
  double rightAlveoliResistanceModifier=1;

  // Currently, There is no way to clear out depositied particulate out of the respiratory systems.
  // Maybe we could have it to cough or some other excretion related method... 
  
  SELiquidSubstanceQuantity* subQ;
  SELiquidSubstanceQuantity* tSubQ;
  const SizeIndependentDepositionEfficencyCoefficient* SIDECoeff;
  double combinedRightBronchodilationEffects = 0.0;
  double combinedLeftBronchodilationEffects = 0.0;
  for (size_t i=0; i<numAerosols; i++)
  {
    //initialize substance
    subQ = m_AerosolMouth->GetSubstanceQuantities()[i];
    //Adjust inflammation coefficient (scaled down):
    inflammationCoefficient = subQ->GetSubstance().GetAerosolization().GetInflammationCoefficient().GetValue();// Once for each subQ
    inflammationCoefficient *= 0.01;
    //Mouth
    SIDECoeff = &m_data.GetSubstances().GetSizeIndependentDepositionEfficencyCoefficient(subQ->GetSubstance());// Once for each subQ
    mouthDepositied_ug = subQ->GetConcentration(MassPerVolumeUnit::ug_Per_mL)*m_AerosolMouth->GetInFlow(VolumePerTimeUnit::mL_Per_s)*m_dt_s*SIDECoeff->GetMouth();
    if (mouthDepositied_ug > subQ->GetMass(MassUnit::ug))
    {
      mouthDepositied_ug = subQ->GetMass(MassUnit::ug);
      subQ->GetMass().SetValue(0, MassUnit::ug);
    }
    else
      subQ->GetMass().IncrementValue(-mouthDepositied_ug, MassUnit::ug); 
    subQ->Balance(BalanceLiquidBy::Mass);
    mouthTotalDepositied_ug = subQ->GetMassDeposited().IncrementValue(mouthDepositied_ug, MassUnit::ug);
    mouthResistanceModifier += mouthTotalDepositied_ug*inflammationCoefficient;
    //Carina
    subQ = m_AerosolCarina->GetSubstanceQuantities()[i];
    carinaDepositied_ug = subQ->GetConcentration(MassPerVolumeUnit::ug_Per_mL)*m_AerosolCarina->GetInFlow(VolumePerTimeUnit::mL_Per_s)*m_dt_s*SIDECoeff->GetCarina();
    if (carinaDepositied_ug > subQ->GetMass(MassUnit::ug))
    {
      carinaDepositied_ug = subQ->GetMass(MassUnit::ug);
      subQ->GetMass().SetValue(0, MassUnit::ug);
    }
    else
      subQ->GetMass().IncrementValue(-carinaDepositied_ug, MassUnit::ug);
    subQ->Balance(BalanceLiquidBy::Mass);
    carinaTotalDepositied_ug = subQ->GetMassDeposited().IncrementValue(carinaDepositied_ug, MassUnit::ug);
    carinaResistanceModifier += carinaTotalDepositied_ug*inflammationCoefficient;
    //Left DeadSpace
    subQ = m_AerosolLeftAnatomicDeadSpace->GetSubstanceQuantities()[i];
    leftDeadSpaceDepositied_ug = subQ->GetConcentration(MassPerVolumeUnit::ug_Per_mL)*m_AerosolLeftAnatomicDeadSpace->GetInFlow(VolumePerTimeUnit::mL_Per_s)*m_dt_s*SIDECoeff->GetDeadSpace();
    if (leftDeadSpaceDepositied_ug > subQ->GetMass(MassUnit::ug))
    {
      leftDeadSpaceDepositied_ug = subQ->GetMass(MassUnit::ug);
      subQ->GetMass().SetValue(0, MassUnit::ug);
    }
    else
      subQ->GetMass().IncrementValue(-leftDeadSpaceDepositied_ug, MassUnit::ug);
    subQ->Balance(BalanceLiquidBy::Mass);
    leftDeadSpaceTotalDepositied_ug = subQ->GetMassDeposited().IncrementValue(leftDeadSpaceDepositied_ug, MassUnit::ug);
    leftDeadSpaceResistanceModifier += leftDeadSpaceTotalDepositied_ug*inflammationCoefficient;
    //Left Alveoli
    subQ = m_AerosolLeftAlveoli->GetSubstanceQuantities()[i];
    leftAlveoliDepositied_ug = subQ->GetConcentration(MassPerVolumeUnit::ug_Per_mL)*m_AerosolLeftAlveoli->GetInFlow(VolumePerTimeUnit::mL_Per_s)*m_dt_s*SIDECoeff->GetAlveoli();
    if (leftAlveoliDepositied_ug > subQ->GetMass(MassUnit::ug))
    {
      leftAlveoliDepositied_ug = subQ->GetMass(MassUnit::ug);
      subQ->GetMass().SetValue(0, MassUnit::ug);
    }
    else
      subQ->GetMass().IncrementValue(-leftAlveoliDepositied_ug, MassUnit::ug);
    subQ->Balance(BalanceLiquidBy::Mass);
    leftAlveoliTotalDepositied_ug = subQ->GetMassDeposited().IncrementValue(leftAlveoliDepositied_ug, MassUnit::ug);
    leftAlveoliResistanceModifier += leftAlveoliTotalDepositied_ug*inflammationCoefficient;
    //Right DeadSpace
    subQ = m_AerosolRightAnatomicDeadSpace->GetSubstanceQuantities()[i];
    rightDeadSpaceDepositied_ug = subQ->GetConcentration(MassPerVolumeUnit::ug_Per_mL)*m_AerosolRightAnatomicDeadSpace->GetInFlow(VolumePerTimeUnit::mL_Per_s)*m_dt_s*SIDECoeff->GetDeadSpace();
    if (rightDeadSpaceDepositied_ug > subQ->GetMass(MassUnit::ug))
    {
      rightDeadSpaceDepositied_ug = subQ->GetMass(MassUnit::ug);
      subQ->GetMass().SetValue(0, MassUnit::ug);
    }
    else
      subQ->GetMass().IncrementValue(-rightDeadSpaceDepositied_ug, MassUnit::ug);
    subQ->Balance(BalanceLiquidBy::Mass);
    rightDeadSpaceTotalDepositied_ug = subQ->GetMassDeposited().IncrementValue(rightDeadSpaceDepositied_ug, MassUnit::ug);
    rightDeadSpaceResistanceModifier += rightDeadSpaceTotalDepositied_ug*inflammationCoefficient;
    //Right Alveoli
    subQ = m_AerosolRightAlveoli->GetSubstanceQuantities()[i];
    rightAlveoliDepositied_ug = subQ->GetConcentration(MassPerVolumeUnit::ug_Per_mL)*m_AerosolRightAlveoli->GetInFlow(VolumePerTimeUnit::mL_Per_s)*m_dt_s*SIDECoeff->GetAlveoli();
    if (rightAlveoliDepositied_ug > subQ->GetMass(MassUnit::ug))
    {
      rightAlveoliDepositied_ug = subQ->GetMass(MassUnit::ug);
      subQ->GetMass().SetValue(0, MassUnit::ug);
    }
    else
      subQ->GetMass().IncrementValue(-rightAlveoliDepositied_ug, MassUnit::ug);
    subQ->Balance(BalanceLiquidBy::Mass);
    rightAlveoliTotalDepositied_ug = subQ->GetMassDeposited().IncrementValue(rightAlveoliDepositied_ug, MassUnit::ug);
    rightAlveoliResistanceModifier += rightAlveoliTotalDepositied_ug*inflammationCoefficient;
    
    // Apply the BronchileModifier dilation effect
    // This is all just tuned to Albuterol - it'll work for other substances, and can be tuned using the other parameters (especially BronchioleModifier)
    if (subQ->GetSubstance().GetState() == eSubstance_State::Liquid)
    {
      // Sum the Bronchiole Effects
      // Must be positive
      double bronchioleModifier = subQ->GetSubstance().GetAerosolization().GetBronchioleModifier().GetValue();

      // Diffuse into Tissues
      // We only process mass deposited on the lungs (dead space and alveoli)
      // We do not currently do anything with the mass in the mouth and carina
      // Could possibly let it go into the stomach somehow... 
      tSubQ = m_LeftLungExtravascular->GetSubstanceQuantity(subQ->GetSubstance());
      tSubQ->GetMass().IncrementValue(leftDeadSpaceDepositied_ug + leftAlveoliDepositied_ug, MassUnit::ug); tSubQ->Balance(BalanceLiquidBy::Mass);    
      combinedLeftBronchodilationEffects += bronchioleModifier * tSubQ->GetConcentration(MassPerVolumeUnit::ug_Per_mL);

      tSubQ = m_RightLungExtravascular->GetSubstanceQuantity(subQ->GetSubstance());
      tSubQ->GetMass().IncrementValue(rightDeadSpaceDepositied_ug + rightAlveoliDepositied_ug, MassUnit::ug); tSubQ->Balance(BalanceLiquidBy::Mass);
      combinedRightBronchodilationEffects += bronchioleModifier * tSubQ->GetConcentration(MassPerVolumeUnit::ug_Per_mL);
    }
  }

  //We're going to make the bronchodilation effects be based off of Albuterol.
  //Experimentally, 1mg of Albuterol given via a spacer device on an endotracheal tube caused a pulmonary resistance change of ~20% @cite mancebo1991effects.
  //The bronchi are ~30% of the total pulmonary resistance, so we'll make a dilation effect really strong with a normal dose.
  //This was tuned using a standard inhaler dose of 90ug
  m_AverageLocalTissueBronchodilationEffects = (combinedLeftBronchodilationEffects + combinedRightBronchodilationEffects) / 2.0;
  combinedLeftBronchodilationEffects = exp(-32894.0 * combinedLeftBronchodilationEffects);
  combinedRightBronchodilationEffects = exp(-32894.0 * combinedRightBronchodilationEffects);  

  // Change resistances due to deposition
  double dTracheaResistance = m_MouthToCarina->GetNextResistance(PressureTimePerVolumeUnit::cmH2O_s_Per_L)*carinaResistanceModifier;
  double dLeftAlveoliResistance = m_LeftAnatomicDeadSpaceToLeftAlveolarDeadSpace->GetNextResistance(PressureTimePerVolumeUnit::cmH2O_s_Per_L)*leftAlveoliResistanceModifier;
  double dRightAlveoliResistance = m_RightAnatomicDeadSpaceToRightAlveolarDeadSpace->GetNextResistance(PressureTimePerVolumeUnit::cmH2O_s_Per_L)*rightAlveoliResistanceModifier;

  double dLeftBronchiResistance = m_CarinaToLeftAnatomicDeadSpace->GetNextResistance(PressureTimePerVolumeUnit::cmH2O_s_Per_L);
  double dRightBronchiResistance = m_CarinaToRightAnatomicDeadSpace->GetNextResistance(PressureTimePerVolumeUnit::cmH2O_s_Per_L);
  dLeftBronchiResistance *= leftDeadSpaceResistanceModifier * combinedLeftBronchodilationEffects;
  dRightBronchiResistance *= rightDeadSpaceResistanceModifier * combinedRightBronchodilationEffects;
  dLeftBronchiResistance = LIMIT(dLeftBronchiResistance, m_RespClosedResistance_cmH2O_s_Per_L, m_RespOpenResistance_cmH2O_s_Per_L);
  dRightBronchiResistance = LIMIT(dRightBronchiResistance, m_RespClosedResistance_cmH2O_s_Per_L, m_RespOpenResistance_cmH2O_s_Per_L);

  m_MouthToCarina->GetNextResistance().SetValue(dTracheaResistance, PressureTimePerVolumeUnit::cmH2O_s_Per_L);
  m_CarinaToLeftAnatomicDeadSpace->GetNextResistance().SetValue(dLeftBronchiResistance, PressureTimePerVolumeUnit::cmH2O_s_Per_L);
  m_CarinaToRightAnatomicDeadSpace->GetNextResistance().SetValue(dRightBronchiResistance, PressureTimePerVolumeUnit::cmH2O_s_Per_L);
  m_LeftAnatomicDeadSpaceToLeftAlveolarDeadSpace->GetNextResistance().SetValue(dLeftAlveoliResistance, PressureTimePerVolumeUnit::cmH2O_s_Per_L);
  m_RightAnatomicDeadSpaceToRightAlveolarDeadSpace->GetNextResistance().SetValue(dRightAlveoliResistance, PressureTimePerVolumeUnit::cmH2O_s_Per_L);
}

//--------------------------------------------------------------------------------------------------
/// \brief
/// Modifies the pressure and/or flow at the mouth
///
/// \details
/// Handles the mechanical ventilation action that adds a flow and pressure source to instantaneously
/// set the respiratory connection (mouth) to user specified values.  
//--------------------------------------------------------------------------------------------------
void Respiratory::MechanicalVentilation()
{
  if (m_data.GetActions().GetPatientActions().HasMechanicalVentilation())
  {
    SEMechanicalVentilation* mv = m_data.GetActions().GetPatientActions().GetMechanicalVentilation();
    // You only get here if action is On
    m_data.SetAirwayMode(eAirwayMode::MechanicalVentilator);

    //Set the substance volume fractions ********************************************
    std::vector<SESubstanceFraction*> gasFractions = mv->GetGasFractions();

    //Reset the substance quantities at the connection
    for (SEGasSubstanceQuantity* subQ : m_MechanicalVentilatorConnection->GetSubstanceQuantities())
      subQ->SetToZero();

    //If no gas fractions specified, assume ambient
    if (gasFractions.empty())
    {
      for (auto s : m_Environment->GetSubstanceQuantities())
      {
        m_MechanicalVentilatorConnection->GetSubstanceQuantity(s->GetSubstance())->GetVolumeFraction().Set(s->GetVolumeFraction());
      }
    }
    else
    {
      //Has fractions defined
      for (auto f : gasFractions)
      {
        SESubstance& sub = f->GetSubstance();
        double fraction = f->GetFractionAmount().GetValue();

        //Do this, just in case it's something new
        m_data.GetSubstances().AddActiveSubstance(sub);

        //Now set it on the connection compartment
        //It has a NaN volume, so this will keep the same volume fraction no matter what's going on around it
        m_MechanicalVentilatorConnection->GetSubstanceQuantity(sub)->GetVolumeFraction().SetValue(fraction);
      }
    }

    //Set the aerosol concentrations ********************************************
    std::vector<SESubstanceConcentration*> liquidConcentrations = mv->GetAerosols();

    //Reset the substance quantities at the connection
    for (SELiquidSubstanceQuantity* subQ : m_MechanicalVentilatorAerosolConnection->GetSubstanceQuantities())
      subQ->SetToZero();

    if (!liquidConcentrations.empty())
    {
      //Has fractions defined
      for (auto f : liquidConcentrations)
      {
        SESubstance& sub = f->GetSubstance();
        SEScalarMassPerVolume concentration = f->GetConcentration();

        //Do this, just in case it's something new
        m_data.GetSubstances().AddActiveSubstance(sub);

        //Now set it on the connection compartment
        //It has a NaN volume, so this will keep the same volume fraction no matter what's going on around it
        m_MechanicalVentilatorAerosolConnection->GetSubstanceQuantity(sub)->GetConcentration().Set(concentration);
      }
    }

    //Apply the instantaneous flow ********************************************
    if (mv->HasFlow())
    {
      m_ConnectionToMouth->GetNextFlowSource().Set(mv->GetFlow());
      //It may or may not be there
      if (!m_ConnectionToMouth->HasFlowSource())
      {
        m_data.GetCircuits().GetRespiratoryAndMechanicalVentilatorCircuit().StateChange();
      }
    }
    else
    {
      //If there's no flow specified, we need to remove the flow source    
      if (m_ConnectionToMouth->HasNextFlowSource())
      {
        m_ConnectionToMouth->GetNextFlowSource().Invalidate();
        m_data.GetCircuits().GetRespiratoryAndMechanicalVentilatorCircuit().StateChange();
      }
    }

    //Apply the instantaneous pressure ********************************************  
    if (mv->HasPressure())
    {
      //This is the pressure above ambient
      m_GroundToConnection->GetNextPressureSource().Set(mv->GetPressure());
    }
    else
    {
      //Pressure is same as ambient
      m_GroundToConnection->GetNextPressureSource().SetValue(0.0, PressureUnit::cmH2O);
    }
  }
  else if (m_data.GetAirwayMode() == eAirwayMode::MechanicalVentilator)
  {
    // Was just turned off
    m_data.SetAirwayMode(eAirwayMode::Free);
    if (m_ConnectionToMouth->HasNextFlowSource())
    {
      m_ConnectionToMouth->GetNextFlowSource().Invalidate();
      m_data.GetCircuits().GetRespiratoryAndMechanicalVentilatorCircuit().StateChange();
    }
  }
}

//--------------------------------------------------------------------------------------------------
/// \brief
/// Applys supplemental oxygen equipment
///
/// \details
/// This will provide supplemental oxygen by connecting and updating the circuit and graph for nasal
/// cannula, simple mask, or nonrebreather mask
//--------------------------------------------------------------------------------------------------
void Respiratory::SupplementalOxygen()
{
  ///\todo - Maybe this and mechanical ventilator should be broken out to their own class, like anesthesia machine?

  if (!m_data.GetActions().GetPatientActions().HasSupplementalOxygen())
    return;
  SESupplementalOxygen* so = m_data.GetActions().GetPatientActions().GetSupplementalOxygen();

  // Get flow
  double flow_L_Per_min = so->GetFlow(VolumePerTimeUnit::L_Per_min);
  //Get tank pressure node and flow control resistor path
  SEFluidCircuitPath* Tank = nullptr;
  SEFluidCircuitPath* OxygenInlet = nullptr;
  SEFluidCircuit* RespirationCircuit = nullptr;
  switch (so->GetDevice())
  {
    case eSupplementalOxygen_Device::None:
    {
      m_data.SetAirwayMode(eAirwayMode::Free);
      m_data.GetActions().GetPatientActions().RemoveSupplementalOxygen();
      return;
    }
    case eSupplementalOxygen_Device::NasalCannula:
    {
      m_data.SetAirwayMode(eAirwayMode::NasalCannula);
      if (!so->HasFlow())
      {
        flow_L_Per_min = 3.5;
        so->GetFlow().SetValue(flow_L_Per_min, VolumePerTimeUnit::L_Per_min);
        Info("Supplemental oxygen flow not set. Using default value of " + std::to_string(flow_L_Per_min) + " L/min.");
      }
      RespirationCircuit = &m_data.GetCircuits().GetActiveRespiratoryCircuit();
      OxygenInlet = RespirationCircuit->GetPath(pulse::NasalCannulaPath::NasalCannulaOxygenInlet);
      Tank = RespirationCircuit->GetPath(pulse::NasalCannulaPath::NasalCannulaPressure);
      break;
    }
    case eSupplementalOxygen_Device::NonRebreatherMask:
    {
      m_data.SetAirwayMode(eAirwayMode::NonRebreatherMask);
      if (!so->HasFlow())
      {
        flow_L_Per_min = 10.0;
        so->GetFlow().SetValue(flow_L_Per_min, VolumePerTimeUnit::L_Per_min);
        Info("Supplemental oxygen flow not set. Using default value of " + std::to_string(flow_L_Per_min) + " L/min.");
      }
      RespirationCircuit = &m_data.GetCircuits().GetActiveRespiratoryCircuit();
      OxygenInlet = RespirationCircuit->GetPath(pulse::NonRebreatherMaskPath::NonRebreatherMaskOxygenInlet);
      Tank = RespirationCircuit->GetPath(pulse::NonRebreatherMaskPath::NonRebreatherMaskPressure);
      break;
    }
    case eSupplementalOxygen_Device::SimpleMask:
    {
      m_data.SetAirwayMode(eAirwayMode::SimpleMask);
      if (!so->HasFlow())
      {
        flow_L_Per_min = 7.5;
        so->GetFlow().SetValue(flow_L_Per_min, VolumePerTimeUnit::L_Per_min);
        Info("Supplemental oxygen flow not set. Using default value of " + std::to_string(flow_L_Per_min) + " L/min.");
      }
      RespirationCircuit = &m_data.GetCircuits().GetActiveRespiratoryCircuit();
      OxygenInlet = RespirationCircuit->GetPath(pulse::SimpleMaskPath::SimpleMaskOxygenInlet);
      Tank = RespirationCircuit->GetPath(pulse::SimpleMaskPath::SimpleMaskPressure);
      break;
    }
    default:
    {
      Error("Ignoring unsupported supplemental oxygen type : " + eSupplementalOxygen_Device_Name(so->GetDevice()));
      return;
    }
  }

  //Use a default tank volume if it wasn't explicity set
  if (!so->HasVolume())
  {
    so->GetVolume().SetValue(425.0, VolumeUnit::L);
    Info("Supplemental oxygen initial tank volume not set. Using default value of 425 L.");
  }
  

  //Decrement volume from the tank
  //Inf volume is assumed to be a wall connection that will never run out
  if (so->GetVolume(VolumeUnit::L) != std::numeric_limits<double>::infinity())
  {
    so->GetVolume().IncrementValue(-flow_L_Per_min * m_dt_min, VolumeUnit::L);
    //Check if the tank is depleated
    if (so->GetVolume(VolumeUnit::L) <= 0.0)
    {
      so->GetVolume().SetValue(0.0, VolumeUnit::L);
      flow_L_Per_min = 0.0;
      /// \event Supplemental Oxygen: Oxygen bottle is exhausted. There is no longer any oxygen to provide.
      m_data.GetEvents().SetEvent(eEvent::SupplementalOxygenBottleExhausted, true, m_data.GetSimulationTime());
    }
    else
      m_data.GetEvents().SetEvent(eEvent::SupplementalOxygenBottleExhausted, false, m_data.GetSimulationTime());
  }
  else
    m_data.GetEvents().SetEvent(eEvent::SupplementalOxygenBottleExhausted, false, m_data.GetSimulationTime());

  //Nonrebreather mask works differently with the bag and doesn't have a pressure source for the tank
  if (so->GetDevice() != eSupplementalOxygen_Device::NonRebreatherMask)
  {
    //Determine flow control resistance
    //Assume pressure outside tank is comparatively approximately ambient
    double tankPressure_cmH2O = Tank->GetNextPressureSource(PressureUnit::cmH2O);
    double resistance_cmH2O_s_Per_L = tankPressure_cmH2O / (flow_L_Per_min / 60.0); //convert from min to s
    resistance_cmH2O_s_Per_L = LIMIT(resistance_cmH2O_s_Per_L, m_DefaultClosedResistance_cmH2O_s_Per_L, m_DefaultOpenResistance_cmH2O_s_Per_L);
    
    //Apply flow
    OxygenInlet->GetNextResistance().SetValue(resistance_cmH2O_s_Per_L, PressureTimePerVolumeUnit::cmH2O_s_Per_L);
  }
  else
  {
    //Handle nonrebreather mask bag
    SEFluidCircuitNode* NonRebreatherMaskBag = RespirationCircuit->GetNode(pulse::NonRebreatherMaskNode::NonRebreatherMaskBag);
    SEFluidCircuitPath* NonRebreatherMaskReservoirValve = RespirationCircuit->GetPath(pulse::NonRebreatherMaskPath::NonRebreatherMaskReservoirValve);

    double flowOut_L_Per_min = 0.0;
    if (NonRebreatherMaskReservoirValve->HasNextFlow())
    {
      flowOut_L_Per_min = NonRebreatherMaskReservoirValve->GetNextFlow(VolumePerTimeUnit::L_Per_min);
    }

    double bagVolume_L = NonRebreatherMaskBag->GetNextVolume(VolumeUnit::L);
    bagVolume_L = bagVolume_L + (flow_L_Per_min - flowOut_L_Per_min) * m_dt_min;
    if (bagVolume_L < 0.0)
    {
      bagVolume_L = 0.0;

      //No air can come from the bag
      OxygenInlet->GetNextResistance().SetValue(m_RespOpenResistance_cmH2O_s_Per_L, PressureTimePerVolumeUnit::cmH2O_s_Per_L);

      /// \event Supplemental Oxygen: The nonrebreather mask is empty. Oxygen may need to be provided at a faster rate.
      m_data.GetEvents().SetEvent(eEvent::NonRebreatherMaskOxygenBagEmpty, true, m_data.GetSimulationTime());
    }
    else
    {
      if (bagVolume_L > 1.0)
      {
        bagVolume_L = 1.0;
      }
      m_data.GetEvents().SetEvent(eEvent::NonRebreatherMaskOxygenBagEmpty, false, m_data.GetSimulationTime());
    }
    
    NonRebreatherMaskBag->GetNextVolume().SetValue(bagVolume_L, VolumeUnit::L);
  }
}

//--------------------------------------------------------------------------------------------------
/// \brief
/// Respiratory driver pressure source
///
/// \details
/// Calculates the muscle pressure source pressure by using the chemical stimuli as feedback control mechanism.
/// The method reads the arterial O2 and CO2 partial pressures. Using these partial pressures, the method calculates the
/// alveolar ventilation from which the method calculates the target breathing frequency. The target breathing frequency is
/// used in the equation for muscle pressure source. The muscle pressure source is used as a driver for ventilation.
/// This method also calculates the drug modifiers that adjusts the depth and frequency of respiration.
//--------------------------------------------------------------------------------------------------
void Respiratory::RespiratoryDriver()
{
  /// \event Patient: Start of exhale/inhale
  if (m_data.GetEvents().IsEventActive(eEvent::StartOfExhale))
    m_data.GetEvents().SetEvent(eEvent::StartOfExhale, false, m_data.GetSimulationTime());
  if (m_data.GetEvents().IsEventActive(eEvent::StartOfInhale))
    m_data.GetEvents().SetEvent(eEvent::StartOfInhale, false, m_data.GetSimulationTime());

  m_BreathingCycleTime_s += m_dt_s;

  //Keep a running average of the Arterial Partial Pressures  
  m_ArterialO2RunningAverage_mmHg->Sample(m_AortaO2->GetPartialPressure(PressureUnit::mmHg));
  m_ArterialCO2RunningAverage_mmHg->Sample(m_AortaCO2->GetPartialPressure(PressureUnit::mmHg));
  //Reset at start of cardiac cycle 
  if (m_data.GetEvents().IsEventActive(eEvent::StartOfCardiacCycle))
  {
    m_ArterialO2PartialPressure_mmHg = m_ArterialO2RunningAverage_mmHg->Value();
    m_ArterialCO2PartialPressure_mmHg = m_ArterialCO2RunningAverage_mmHg->Value();

    m_ArterialO2RunningAverage_mmHg->Clear();
    m_ArterialCO2RunningAverage_mmHg->Clear();
  }

#ifdef TUNING
  m_ArterialO2PartialPressure_mmHg = 95.0;
  m_ArterialCO2PartialPressure_mmHg = 40.0;
#endif  

  if (m_ConsciousBreathing) //Conscious breathing 
  {
    SEConsciousRespiration* cr = m_data.GetActions().GetPatientActions().GetConsciousRespiration();
    SEConsciousRespirationCommand* cmd = cr->GetActiveCommand();
    SEUseInhaler* ui = dynamic_cast<SEUseInhaler*>(cmd);
    if (ui != nullptr)
    {
      //We're using the inhaler, so just wait at this driver pressure
      m_DriverPressure_cmH2O = m_DriverPressurePath->GetPressureSource(PressureUnit::cmH2O);
      m_ConsciousBreathing = false;
    }
    else
    {
      //Just increase/decrease the driver pressure linearly to achieve the desired 
      //pressure (attempting to reach a specific volume) at the end of the user specified period
      double Time_s = m_ConsciousRespirationPeriod_s - m_ConsciousRespirationRemainingPeriod_s;
      double Slope = (m_ConsciousEndPressure_cmH2O - m_ConsciousStartPressure_cmH2O) / m_ConsciousRespirationPeriod_s;
      //Form of y=mx+b
      m_DriverPressure_cmH2O = Slope * Time_s + m_ConsciousStartPressure_cmH2O;

      //Make it so we start a fresh cycle when we go back to spontaneous breathing
      //Adding 2.0 * timestamp just makes it greater than the TotalBreathingCycleTime_s
      m_VentilationFrequency_Per_min = m_data.GetCurrentPatient().GetRespirationRateBaseline(FrequencyUnit::Per_min);
      m_BreathingCycleTime_s = 60.0 / m_VentilationFrequency_Per_min + 2.0 * m_dt_s;
    }
  }
  else //Spontaneous (i.e. not conscious) breathing
  {
    double TotalBreathingCycleTime_s = 0.0;
    if (m_VentilationFrequency_Per_min < 1.0)
    {
      TotalBreathingCycleTime_s = 0.0;
    }
    else
    {
      TotalBreathingCycleTime_s = 60.0 / m_VentilationFrequency_Per_min; //Total time of one breathing cycle  
    }

    //Prepare for the next cycle -------------------------------------------------------------------------------
    if (m_BreathingCycleTime_s > TotalBreathingCycleTime_s) //End of the cycle or currently not breathing
    {
      // Make a cardicArrestEffect that is 1.0 unless cardiac arrest is true
      double cardiacArrestEffect = 1.0;
      // If the cv system parameter is true, then make the cardicArrestEffect = 0
      if (m_data.GetEvents().IsEventActive(eEvent::CardiacArrest))
      {
        cardiacArrestEffect = 0.0;
      }

      //Ventilatory Negative Feedback Control *************************************************************************
      double PeripheralCO2PartialPressureSetPoint = m_data.GetConfiguration().GetPeripheralControllerCO2PressureSetPoint(PressureUnit::mmHg);
      double CentralCO2PartialPressureSetPoint = m_data.GetConfiguration().GetCentralControllerCO2PressureSetPoint(PressureUnit::mmHg);

      double metabolicModifier = 1.0;
      double TMR_W = m_data.GetEnergy().GetTotalMetabolicRate(PowerUnit::W);
      double BMR_W = m_data.GetCurrentPatient().GetBasalMetabolicRate(PowerUnit::W);
      double metabolicFraction = TMR_W / BMR_W;

      //Get Drug Effects
      SEDrugSystem& Drugs = m_data.GetDrugs();
      double DrugRRChange_Per_min = Drugs.GetRespirationRateChange(FrequencyUnit::Per_min);
      double DrugsTVChange_L = Drugs.GetTidalVolumeChange(VolumeUnit::L);
      double NMBModifier = 1.0;
      double SedationModifier = 1.0;
      //Make changes to Respiration Rate change based on the neuromuscular block level
      if (Drugs.GetNeuromuscularBlockLevel().GetValue() > 0.135)
      {
        NMBModifier = 0.0;
        DrugRRChange_Per_min = 0.0;
        DrugsTVChange_L = 0.0;
      }
      else if (Drugs.GetNeuromuscularBlockLevel().GetValue() > 0.11)
      {
        NMBModifier = 0.5;
        DrugRRChange_Per_min = 0.0;
        DrugsTVChange_L = 0.0;
      }
      else if (Drugs.GetSedationLevel().GetValue() > 0.15 && std::abs(m_data.GetCurrentPatient().GetRespirationRateBaseline(FrequencyUnit::Per_min) + DrugRRChange_Per_min) < 1.0)
      {
        SedationModifier = 0.0;
        DrugRRChange_Per_min = 0.0;
        DrugsTVChange_L = 0.0;
      }

      //Calculate the target Alveolar Ventilation based on the Arterial O2 and CO2 concentrations
      double dTargetAlveolarVentilation_L_Per_min = m_PeripheralControlGainConstant * exp(-0.05*m_ArterialO2PartialPressure_mmHg)*MAX(0., m_ArterialCO2PartialPressure_mmHg - PeripheralCO2PartialPressureSetPoint); //Peripheral portion
      dTargetAlveolarVentilation_L_Per_min += m_CentralControlGainConstant * MAX(0., m_ArterialCO2PartialPressure_mmHg - CentralCO2PartialPressureSetPoint); //Central portion

      //Metabolic modifier is used to drive the system to reasonable levels achievable during increased metabolic exertion
      //The modifier is tuned to achieve the correct respiratory response for near maximal exercise. A linear relationship is assumed
      // for the respiratory effects due to increased metabolic exertion    
      double tunedVolumeMetabolicSlope = 0.2; //Tuned fractional increase of the tidal volume due to increased metabolic rate
      metabolicModifier = 1.0 + tunedVolumeMetabolicSlope * (metabolicFraction - 1.0);
      dTargetAlveolarVentilation_L_Per_min *= metabolicModifier;

      //Only move it part way to where it wants to be.
      //If it stays there, it will get there, just more controlled/slowly.
      //This is needed so we don't overshoot and introduce low frequency oscillations into the system (i.e. overdamped).
      double targetCO2PP_mmHg = 40.0;
      double changeFraction = 1.0;
      if (m_data.GetState() <= EngineState::InitialStabilization)
      {
        //This gets it nice and stable
        changeFraction = 0.1;
      }
      else
      {
        //This keeps it from going crazy with modifiers applied
        changeFraction = std::abs(m_ArterialCO2PartialPressure_mmHg - targetCO2PP_mmHg) / targetCO2PP_mmHg * 0.5;
      }
      changeFraction = MIN(changeFraction, 1.0);

#ifdef TUNING
      changeFraction = 1.0;
#endif

      dTargetAlveolarVentilation_L_Per_min = m_PreviousTargetAlveolarVentilation_L_Per_min + (dTargetAlveolarVentilation_L_Per_min - m_PreviousTargetAlveolarVentilation_L_Per_min) * changeFraction;
      m_PreviousTargetAlveolarVentilation_L_Per_min = dTargetAlveolarVentilation_L_Per_min;

      //Target Tidal Volume (i.e. Driver amplitude) *************************************************************************
      //Calculate the target Tidal Volume based on the Alveolar Ventilation
      double dTargetPulmonaryVentilation_L_Per_min = dTargetAlveolarVentilation_L_Per_min;

      double dMaximumPulmonaryVentilationRate = m_data.GetConfiguration().GetPulmonaryVentilationRateMaximum(VolumePerTimeUnit::L_Per_min);
      /// \event Patient: Maximum Pulmonary Ventilation Rate : Pulmonary ventilation exceeds maximum value
      if (dTargetPulmonaryVentilation_L_Per_min > dMaximumPulmonaryVentilationRate)
      {
        dTargetPulmonaryVentilation_L_Per_min = dMaximumPulmonaryVentilationRate;
        m_data.GetEvents().SetEvent(eEvent::MaximumPulmonaryVentilationRate, true, m_data.GetSimulationTime());
      }

      if (dTargetPulmonaryVentilation_L_Per_min < dMaximumPulmonaryVentilationRate && m_data.GetEvents().IsEventActive(eEvent::MaximumPulmonaryVentilationRate))
      {
        m_data.GetEvents().SetEvent(eEvent::MaximumPulmonaryVentilationRate, false, m_data.GetSimulationTime());
      }

      //Calculate the target Tidal Volume based on the calculated target pulmonary ventilation, plot slope (determined during initialization), and x-intercept
      double dTargetTidalVolume_L = dTargetPulmonaryVentilation_L_Per_min / m_VentilationToTidalVolumeSlope + m_VentilationTidalVolumeIntercept;

      //Modify the target tidal volume due to other external effects - probably eventually replaced by the Nervous system
      dTargetTidalVolume_L *= cardiacArrestEffect * NMBModifier;

      //Apply Drug Effects to the target tidal volume
      dTargetTidalVolume_L += DrugsTVChange_L;

      //This is a piecewise function that plateaus at the Tidal Volume equal to 1/2 * Vital Capacity
      //The Respiration Rate will make up for the Alveoli Ventilation difference
      double dHalfVitalCapacity_L = m_data.GetCurrentPatient().GetVitalCapacity(VolumeUnit::L) / 2;
      dTargetTidalVolume_L = MIN(dTargetTidalVolume_L, dHalfVitalCapacity_L);

      //Map the Target Tidal Volume to the Driver
      double TargetVolume_L = m_data.GetCurrentPatient().GetFunctionalResidualCapacity(VolumeUnit::L) + dTargetTidalVolume_L;
      m_PeakRespiratoryDrivePressure_cmH2O = VolumeToDriverPressure(TargetVolume_L);
      //There's a maximum force the driver can try to achieve
      m_PeakRespiratoryDrivePressure_cmH2O = MAX(m_PeakRespiratoryDrivePressure_cmH2O, m_MaxDriverPressure_cmH2O);


      //Respiration Rate (i.e. Driver frequency) *************************************************************************
      //Calculate the Respiration Rate given the Alveolar Ventilation and the Target Tidal Volume
      if (dTargetTidalVolume_L == 0.0) //Can't divide by zero
      {
        m_VentilationFrequency_Per_min = 0.0;
        m_NotBreathing = true;
      }
      else
      {
        m_VentilationFrequency_Per_min = dTargetPulmonaryVentilation_L_Per_min / dTargetTidalVolume_L; //breaths/min
        m_VentilationFrequency_Per_min *= NMBModifier * SedationModifier;
        m_VentilationFrequency_Per_min += DrugRRChange_Per_min;
        m_NotBreathing = false;
      }

      m_VentilationFrequency_Per_min = LIMIT(m_VentilationFrequency_Per_min, 0.0, dMaximumPulmonaryVentilationRate / dHalfVitalCapacity_L);

      //Patient Definition *************************************************************************
      //We need to hit the patient's defined Respiration Rate Baseline, no matter what,
      //so we'll keep adjusting the slope of the function to achieve this
      if (m_data.GetState() <= EngineState::InitialStabilization)
      {
        double dRespirationRateBaseline_Per_min = m_data.GetCurrentPatient().GetRespirationRateBaseline(FrequencyUnit::Per_min);
        double dPercentError = (m_VentilationFrequency_Per_min - dRespirationRateBaseline_Per_min) / dRespirationRateBaseline_Per_min; //negative if too low

        //Amplitude set-point - this will set the Tidal Volume baseline when O2 and CO2 are at the correct/balanced level
        m_VentilationToTidalVolumeSlope = m_VentilationToTidalVolumeSlope - m_VentilationToTidalVolumeSlope * dPercentError;

        //Put bounds on this
        m_VentilationToTidalVolumeSlope = LIMIT(m_VentilationToTidalVolumeSlope, 1.0, 100.0);
      }

#ifdef TUNING
      m_VentilationToTidalVolumeSlope = 40.0;
      m_VentilationFrequency_Per_min = 16.0;
#endif

      SetBreathCycleFractions();

      m_BreathingCycleTime_s = 0.0;

      //KEEP COMMENTED OUT - this is for keeping the driver constant while debugging
      //m_VentilationFrequency_Per_min = 16.0;
      //m_PeakRespiratoryDrivePressure = -11.3;

      m_data.GetEvents().SetEvent(eEvent::StartOfInhale, true, m_data.GetSimulationTime()); ///\todo rename "StartOfConsiousInhale"
    }

    //Run the driver based on the waveform -------------------------------------------------------------------------------
    TotalBreathingCycleTime_s = 60.0 / m_VentilationFrequency_Per_min; //Recalculate, in case we just calculated a new breath

    double PeakExpiratoryPressure = 0.0; //For potential future implementation

    double InspiratoryRiseTimeStart_s = 0.0;
    double InspiratoryHoldTimeStart_s = InspiratoryRiseTimeStart_s + m_InspiratoryRiseFraction * TotalBreathingCycleTime_s;
    double InspiratoryReleaseTimeStart_s = InspiratoryHoldTimeStart_s + m_InspiratoryHoldFraction * TotalBreathingCycleTime_s;
    double InspiratoryToExpiratoryPauseTimeStart_s = InspiratoryReleaseTimeStart_s + m_InspiratoryReleaseFraction * TotalBreathingCycleTime_s;
    double ExpiratoryRiseTimeStart_s = InspiratoryToExpiratoryPauseTimeStart_s + m_InspiratoryToExpiratoryPauseFraction * TotalBreathingCycleTime_s;
    double ExpiratoryHoldTimeStart_s = ExpiratoryRiseTimeStart_s + m_ExpiratoryRiseFraction * TotalBreathingCycleTime_s;
    double ExpiratoryReleaseTimeStart_s = ExpiratoryHoldTimeStart_s + m_ExpiratoryHoldFraction * TotalBreathingCycleTime_s;
    double ResidueFractionTimeStart_s = ExpiratoryReleaseTimeStart_s + m_ExpiratoryReleaseFraction * TotalBreathingCycleTime_s;

    m_DriverInspirationTime_s = ExpiratoryRiseTimeStart_s;


    if (m_BreathingCycleTime_s >= InspiratoryReleaseTimeStart_s &&
      m_BreathingCycleTime_s < InspiratoryReleaseTimeStart_s + m_dt_s) //Only call this once per cycle