Respiratory.cpp 158 KB
Newer Older
Aaron Bray's avatar
Aaron Bray committed
1
2
/* Distributed under the Apache License, Version 2.0.
   See accompanying NOTICE file for details.*/
Aaron Bray's avatar
Aaron Bray committed
3
4

#include "stdafx.h"
5
6
7
8
#include "physiology/Respiratory.h"
#include "controller/Circuits.h"
#include "controller/Compartments.h"
#include "controller/Substances.h"
9
10
#include "PulseConfiguration.h"
// Conditions 
11
#include "engine/SEConditionManager.h"
12
13
14
#include "patient/conditions/SEChronicObstructivePulmonaryDisease.h"
#include "patient/conditions/SELobarPneumonia.h"
#include "patient/conditions/SEImpairedAlveolarExchange.h"
Aaron Bray's avatar
Aaron Bray committed
15
#include "patient/conditions/SEAcuteRespiratoryDistressSyndrome.h"
16
// Actions
17
18
#include "engine/SEActionManager.h"
#include "engine/SEPatientActionCollection.h"
Aaron Bray's avatar
Aaron Bray committed
19
#include "patient/actions/SEAcuteRespiratoryDistressSyndromeExacerbation.h"
20
21
22
23
#include "patient/actions/SEAirwayObstruction.h"
#include "patient/actions/SEAsthmaAttack.h"
#include "patient/actions/SEBronchoconstriction.h"
#include "patient/actions/SEChestOcclusiveDressing.h"
Aaron Bray's avatar
Aaron Bray committed
24
#include "patient/actions/SEChronicObstructivePulmonaryDiseaseExacerbation.h"
25
26
27
28
29
#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"
30
#include "patient/actions/SEDyspnea.h"
31
#include "patient/actions/SEIntubation.h"
Aaron Bray's avatar
Aaron Bray committed
32
#include "patient/actions/SELobarPneumoniaExacerbation.h"
33
34
#include "patient/actions/SEMechanicalVentilation.h"
#include "patient/actions/SENeedleDecompression.h"
35
#include "patient/actions/SERespiratoryFatigue.h"
36
#include "patient/actions/SESupplementalOxygen.h"
37
38
39
40
41
42
43
44
45
#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"
46
#include "engine/SEEventManager.h"
47
48
#include "substance/SESubstance.h"
#include "substance/SESubstanceAerosolization.h"
49
#include "substance/SESubstanceConcentration.h"
50
51
52
53
#include "substance/SESubstanceFraction.h"
#include "substance/SESubstanceTransport.h"
#include "circuit/fluid/SEFluidCircuitCalculator.h"
#include "circuit/fluid/SEFluidCircuit.h"
Aaron Bray's avatar
Aaron Bray committed
54
#include "circuit/fluid/SEFluidCircuit.h"
55
56
57
58
#include "circuit/fluid/SEFluidCircuitNode.h"
#include "circuit/fluid/SEFluidCircuitPath.h"
#include "compartment/fluid/SEGasCompartmentGraph.h"
#include "compartment/fluid/SELiquidCompartmentGraph.h"
Aaron Bray's avatar
Aaron Bray committed
59
60
// Properties
#include "properties/SEScalar0To1.h"
61
62
#include "properties/SEScalarArea.h"
#include "properties/SEScalarFrequency.h"
Aaron Bray's avatar
Aaron Bray committed
63
#include "properties/SEScalarInversePressure.h"
64
#include "properties/SEScalarInverseVolume.h"
Aaron Bray's avatar
Aaron Bray committed
65
66
67
#include "properties/SEScalarLength.h"
#include "properties/SEScalarMass.h"
#include "properties/SEScalarMassPerVolume.h"
68
69
70
#include "properties/SEScalarNegative1To1.h"
#include "properties/SEScalarPower.h"
#include "properties/SEScalarPressure.h"
71
#include "properties/SEScalarPressurePerVolume.h"
72
#include "properties/SEScalarPressureTimePerVolume.h"
73
74
75
#include "properties/SEScalarTime.h"
#include "properties/SEScalarVolume.h"
#include "properties/SEScalarVolumePerTime.h"
76
#include "properties/SEScalarVolumePerPressure.h"
Aaron Bray's avatar
Aaron Bray committed
77
#include "properties/SEFunctionVolumeVsTime.h"
78
#include "properties/SERunningAverage.h"
79
// Data Track
80
81
#include "utils/DataTrack.h"

Aaron Bray's avatar
Aaron Bray committed
82
83
84
85
#ifdef _MSC_VER 
#pragma warning( disable : 4305 4244 )  // Disable some warning messages
#endif

86
87
88
89
//Flag for data tracks
//Should be commented out, unless debugging
//#define DEBUG

Aaron Bray's avatar
Aaron Bray committed
90
91
//Flag for setting things constant to test
//Should be commented out, unless debugging/tuning
Jeff Webb's avatar
Jeff Webb committed
92
//#define TUNING
Aaron Bray's avatar
Aaron Bray committed
93

94
Respiratory::Respiratory(PulseController& data) : SERespiratorySystem(data.GetLogger()), m_data(data)
Aaron Bray's avatar
Aaron Bray committed
95
{
96
97
98
  m_BloodPHRunningAverage = new SERunningAverage();
  m_ArterialO2RunningAverage_mmHg = new SERunningAverage();
  m_ArterialCO2RunningAverage_mmHg = new SERunningAverage();
99

100
  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());
101
102
  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());
Aaron Bray's avatar
Aaron Bray committed
103
  Clear();
Aaron Bray's avatar
Aaron Bray committed
104
105
106
107
108
}

Respiratory::~Respiratory()
{
  Clear();
109
110
111
112
113
114
115
  delete m_BloodPHRunningAverage;
  delete m_ArterialO2RunningAverage_mmHg;
  delete m_ArterialCO2RunningAverage_mmHg;

  delete m_Calculator;
  delete m_GasTransporter;
  delete m_AerosolTransporter;
Aaron Bray's avatar
Aaron Bray committed
116
117
118
119
120
121
122
123
124
125
}

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

  m_Environment = nullptr;
  m_AerosolMouth = nullptr;
  m_AerosolCarina = nullptr;
126
127
  m_AerosolLeftAnatomicDeadSpace = nullptr;
  m_AerosolLeftAlveolarDeadSpace = nullptr;
Aaron Bray's avatar
Aaron Bray committed
128
  m_AerosolLeftAlveoli = nullptr;
129
130
  m_AerosolRightAnatomicDeadSpace = nullptr;
  m_AerosolRightAlveolarDeadSpace = nullptr;
Aaron Bray's avatar
Aaron Bray committed
131
132
133
134
135
136
137
138
  m_AerosolRightAlveoli = nullptr;
  m_Lungs = nullptr;
  m_LeftLungExtravascular = nullptr;
  m_RightLungExtravascular = nullptr;
  m_Carina = nullptr;
  m_AortaO2 = nullptr;
  m_AortaCO2 = nullptr;
  m_MechanicalVentilatorConnection = nullptr;
139
  m_MechanicalVentilatorAerosolConnection = nullptr;
140
  m_PleuralCavity = nullptr;
Aaron Bray's avatar
Aaron Bray committed
141
142
143

  m_RespiratoryCircuit = nullptr;

Aaron Bray's avatar
Aaron Bray committed
144
  m_Mouth = nullptr;
Aaron Bray's avatar
Aaron Bray committed
145
  m_LeftAlveoli = nullptr;
146
147
  m_LeftAnatomicDeadSpace = nullptr;
  m_LeftAlveolarDeadSpace = nullptr;
Aaron Bray's avatar
Aaron Bray committed
148
149
150
  m_LeftPleural = nullptr;
  m_RespiratoryMuscle = nullptr;
  m_RightAlveoli = nullptr;
151
152
  m_RightAnatomicDeadSpace = nullptr;
  m_RightAlveolarDeadSpace = nullptr;
Aaron Bray's avatar
Aaron Bray committed
153
154
155
156
  m_RightPleural = nullptr;

  m_CarinaToLeftAnatomicDeadSpace = nullptr;
  m_CarinaToRightAnatomicDeadSpace = nullptr;
157
158
159
160
  m_LeftAnatomicDeadSpaceToLeftAlveolarDeadSpace = nullptr;
  m_RightAnatomicDeadSpaceToRightAlveolarDeadSpace = nullptr;
  m_LeftAlveolarDeadSpaceToLeftAlveoli = nullptr;
  m_RightAlveolarDeadSpaceToRightAlveoli = nullptr;
Aaron Bray's avatar
Aaron Bray committed
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
  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;

179
180
181
  m_BloodPHRunningAverage->Clear();
  m_ArterialO2RunningAverage_mmHg->Clear();
  m_ArterialCO2RunningAverage_mmHg->Clear();
Aaron Bray's avatar
Aaron Bray committed
182
183
184
185
186
187
188
189
}

//--------------------------------------------------------------------------------------------------
/// \brief
/// Initializes system properties to valid homeostatic values.
//--------------------------------------------------------------------------------------------------
void Respiratory::Initialize()
{
190
  PulseSystem::Initialize();
Aaron Bray's avatar
Aaron Bray committed
191
192

  //Vital signs
193
  m_NotBreathing = false;
Aaron Bray's avatar
Aaron Bray committed
194
195
  m_TopBreathTotalVolume_L = 0.0;
  m_TopBreathAlveoliVolume_L = 0.0;
196
  m_TopBreathPleuralVolume_L = 0.0;
Aaron Bray's avatar
Aaron Bray committed
197
  m_TopBreathPleuralPressure_cmH2O = 0.0;
198
199
  m_TopBreathAlveoliPressure_cmH2O = 0.0;
  m_TopBreathDriverPressure_cmH2O = 0.0;
Aaron Bray's avatar
Aaron Bray committed
200
  m_LastCardiacCycleBloodPH = 7.4;
201
202
203
204
205
  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;
206
  m_BottomBreathPleuralVolume_L = 0.0;
207
  m_BottomBreathPleuralPressure_cmH2O = 0.0;
208
209
210
211
  m_BottomBreathAlveoliPressure_cmH2O = 0.0;
  m_BottomBreathDriverPressure_cmH2O = 0.0;
  m_PeakAlveolarPressure_cmH2O = 0.0;
  m_MaximalAlveolarPressure_cmH2O = 0.0;
Aaron Bray's avatar
Aaron Bray committed
212
213

  //Driver
214
  m_MaxDriverPressure_cmH2O = -50.0;
Aaron Bray's avatar
Aaron Bray committed
215
  m_ElapsedBreathingCycleTime_min = 0.0;
216
  m_TopBreathElapsedTime_min = 0.0;
Aaron Bray's avatar
Aaron Bray committed
217
218
  m_BreathingCycle = false;
  m_BreathingCycleTime_s = 0.0;
219
  m_VentilationFrequency_Per_min = m_data.GetCurrentPatient().GetRespirationRateBaseline(FrequencyUnit::Per_min);
220
221
  m_DriverPressure_cmH2O = 0.0;
  m_DriverPressureMin_cmH2O = 0.0;
Aaron Bray's avatar
Aaron Bray committed
222
223
  m_VentilationToTidalVolumeSlope = 30.0;
  //The peak driver pressure is the pressure above the default pressure
224
  m_PeakRespiratoryDrivePressure_cmH2O = 0.0;
Aaron Bray's avatar
Aaron Bray committed
225
226
  m_ArterialO2PartialPressure_mmHg = m_AortaO2->GetPartialPressure(PressureUnit::mmHg);
  m_ArterialCO2PartialPressure_mmHg = m_AortaCO2->GetPartialPressure(PressureUnit::mmHg);
227
  m_PreviousTargetAlveolarVentilation_L_Per_min = m_data.GetCurrentPatient().GetTidalVolumeBaseline(VolumeUnit::L) * m_VentilationFrequency_Per_min;
Aaron Bray's avatar
Aaron Bray committed
228
229
  m_AverageLocalTissueBronchodilationEffects = 0.0;

230
  m_IERatioScaleFactor = 1.0;
Aaron Bray's avatar
Aaron Bray committed
231
232
233
234
235
236
237
238
239
240
241

  // 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
242
243
  double TidalVolume_L = m_data.GetCurrentPatient().GetTidalVolumeBaseline(VolumeUnit::L);
  double RespirationRate_Per_min = m_data.GetCurrentPatient().GetRespirationRateBaseline(FrequencyUnit::Per_min);
Aaron Bray's avatar
Aaron Bray committed
244
245
  GetTidalVolume().SetValue(TidalVolume_L, VolumeUnit::L);
  GetRespirationRate().SetValue(RespirationRate_Per_min, FrequencyUnit::Per_min);
246
  GetInspiratoryExpiratoryRatio().SetValue(0.5);
Aaron Bray's avatar
Aaron Bray committed
247
248
  GetCarricoIndex().SetValue(452.0, PressureUnit::mmHg);

249
250
  double AnatomicDeadSpace_L = m_LeftAnatomicDeadSpace->GetVolumeBaseline(VolumeUnit::L) + m_RightAnatomicDeadSpace->GetVolumeBaseline(VolumeUnit::L);
  double AlveolarDeadSpace_L = 0.0;
251
252
  GetAnatomicDeadSpace().SetValue(AnatomicDeadSpace_L, VolumeUnit::L);
  GetAlveolarDeadSpace().SetValue(AlveolarDeadSpace_L, VolumeUnit::L);
253
  GetPhysiologicDeadSpace().SetValue(AnatomicDeadSpace_L + AlveolarDeadSpace_L, VolumeUnit::L);
254
255

  GetTotalAlveolarVentilation().SetValue(RespirationRate_Per_min * TidalVolume_L, VolumePerTimeUnit::L_Per_min);
Aaron Bray's avatar
Aaron Bray committed
256
  GetTotalPulmonaryVentilation().SetValue(RespirationRate_Per_min * TidalVolume_L, VolumePerTimeUnit::L_Per_min);
257
  GetTotalDeadSpaceVentilation().SetValue(AnatomicDeadSpace_L * RespirationRate_Per_min, VolumePerTimeUnit::L_Per_min);
Aaron Bray's avatar
Aaron Bray committed
258
  GetSpecificVentilation().SetValue(0.21);
259

Aaron Bray's avatar
Aaron Bray committed
260
261
  GetEndTidalCarbonDioxideFraction().SetValue(0.0827);
  GetEndTidalCarbonDioxidePressure().SetValue(0.0, PressureUnit::mmHg);
262
263
264
265
266
267
268
269
270
271
272
  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);

273
  GetPulmonaryCompliance().SetValue(0.1, VolumePerPressureUnit::L_Per_cmH2O);
274
  GetLungCompliance().SetValue(0.1, VolumePerPressureUnit::L_Per_cmH2O);
275
  GetChestWallCompliance().SetValue(0.2, VolumePerPressureUnit::L_Per_cmH2O);
276
  GetPulmonaryElastance().SetValue(1.0 / 0.1, PressurePerVolumeUnit::cmH2O_Per_L);
Aaron Bray's avatar
Aaron Bray committed
277

278
279
280
281
282
283
284
285
286
287
  // 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;

Aaron Bray's avatar
Aaron Bray committed
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
  //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
307
308
309
310
  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);
Aaron Bray's avatar
Aaron Bray committed
311
312
313
314
315
  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);
316
317
  m_MinimumAllowableTidalVolume_L = m_data.GetConfiguration().GetMinimumAllowableTidalVolume(VolumeUnit::L);
  m_MinimumAllowableInpiratoryAndExpiratoryPeriod_s = m_data.GetConfiguration().GetMinimumAllowableInpiratoryAndExpiratoryPeriod(TimeUnit::s);
Aaron Bray's avatar
Aaron Bray committed
318
  //Compartments
Aaron Bray's avatar
Aaron Bray committed
319
320
321
  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);
322
323
  m_AerosolLeftAnatomicDeadSpace = m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::LeftAnatomicDeadSpace);
  m_AerosolLeftAlveolarDeadSpace = m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::LeftAlveolarDeadSpace);
Aaron Bray's avatar
Aaron Bray committed
324
  m_AerosolLeftAlveoli = m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::LeftAlveoli);
325
326
  m_AerosolRightAnatomicDeadSpace = m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::RightAnatomicDeadSpace);
  m_AerosolRightAlveolarDeadSpace = m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::RightAlveolarDeadSpace);
Aaron Bray's avatar
Aaron Bray committed
327
328
  m_AerosolRightAlveoli = m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::RightAlveoli);
  m_Lungs = m_data.GetCompartments().GetGasCompartment(pulse::PulmonaryCompartment::Lungs);
329
330
  m_LeftLung = m_data.GetCompartments().GetGasCompartment(pulse::PulmonaryCompartment::LeftLung);
  m_RightLung = m_data.GetCompartments().GetGasCompartment(pulse::PulmonaryCompartment::RightLung);
331
  m_PleuralCavity = m_data.GetCompartments().GetGasCompartment(pulse::PulmonaryCompartment::PleuralCavity);
Aaron Bray's avatar
Aaron Bray committed
332
333
334
  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);
Aaron Bray's avatar
Aaron Bray committed
335
  m_CarinaO2 = m_Carina->GetSubstanceQuantity(m_data.GetSubstances().GetO2());
Aaron Bray's avatar
Aaron Bray committed
336
  SELiquidCompartment* Aorta = m_data.GetCompartments().GetLiquidCompartment(pulse::VascularCompartment::Aorta);
Aaron Bray's avatar
Aaron Bray committed
337
338
  m_AortaO2 = Aorta->GetSubstanceQuantity(m_data.GetSubstances().GetO2());
  m_AortaCO2 = Aorta->GetSubstanceQuantity(m_data.GetSubstances().GetCO2());
Aaron Bray's avatar
Aaron Bray committed
339
340
341
  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);
342
  m_MechanicalVentilatorAerosolConnection = m_data.GetCompartments().GetLiquidCompartment(pulse::MechanicalVentilatorCompartment::Connection);
Aaron Bray's avatar
Aaron Bray committed
343
  // Compartments we will process aerosol effects on
Aaron Bray's avatar
Aaron Bray committed
344
  m_AerosolEffects.clear();
Aaron Bray's avatar
Aaron Bray committed
345
346
  m_AerosolEffects.push_back(m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::Carina));
  m_AerosolEffects.push_back(m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::LeftAlveoli));
347
348
  m_AerosolEffects.push_back(m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::LeftAnatomicDeadSpace));
  m_AerosolEffects.push_back(m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::LeftAlveolarDeadSpace));
Aaron Bray's avatar
Aaron Bray committed
349
  m_AerosolEffects.push_back(m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::RightAlveoli));
350
351
  m_AerosolEffects.push_back(m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::RightAnatomicDeadSpace));
  m_AerosolEffects.push_back(m_data.GetCompartments().GetLiquidCompartment(pulse::PulmonaryCompartment::RightAlveolarDeadSpace));
Aaron Bray's avatar
Aaron Bray committed
352
353
354
  //Circuits
  m_RespiratoryCircuit = &m_data.GetCircuits().GetRespiratoryCircuit();
  //Nodes
355
  m_Mouth = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::Mouth);
Aaron Bray's avatar
Aaron Bray committed
356
  m_LeftAlveoli = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::LeftAlveoli);
357
358
  m_LeftAnatomicDeadSpace = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::LeftAnatomicDeadSpace);
  m_LeftAlveolarDeadSpace = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::LeftAlveolarDeadSpace);
Aaron Bray's avatar
Aaron Bray committed
359
360
361
  m_LeftPleural = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::LeftPleural);
  m_RespiratoryMuscle = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::RespiratoryMuscle);
  m_RightAlveoli = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::RightAlveoli);
362
363
  m_RightAnatomicDeadSpace = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::RightAnatomicDeadSpace);
  m_RightAlveolarDeadSpace = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::RightAlveolarDeadSpace);
Aaron Bray's avatar
Aaron Bray committed
364
365
366
  m_RightPleural = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::RightPleural);
  m_Ambient = m_RespiratoryCircuit->GetNode(pulse::EnvironmentNode::Ambient);
  m_Stomach = m_RespiratoryCircuit->GetNode(pulse::RespiratoryNode::Stomach);
Aaron Bray's avatar
Aaron Bray committed
367
  //Paths
Aaron Bray's avatar
Aaron Bray committed
368
369
  m_CarinaToLeftAnatomicDeadSpace = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::CarinaToLeftAnatomicDeadSpace);
  m_CarinaToRightAnatomicDeadSpace = m_RespiratoryCircuit->GetPath(pulse::RespiratoryPath::CarinaToRightAnatomicDeadSpace);
370
371
372
373
  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);
Aaron Bray's avatar
Aaron Bray committed
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
  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);
Aaron Bray's avatar
Aaron Bray committed
389
390

  /// \todo figure out how to modify these resistances without getting the cv circuit - maybe add a parameter, like baroreceptors does
Aaron Bray's avatar
Aaron Bray committed
391
392
  m_RightPulmonaryCapillary = m_data.GetCircuits().GetCardiovascularCircuit().GetPath(pulse::CardiovascularPath::RightPulmonaryCapillariesToRightPulmonaryVeins);
  m_LeftPulmonaryCapillary = m_data.GetCircuits().GetCardiovascularCircuit().GetPath(pulse::CardiovascularPath::LeftPulmonaryCapillariesToLeftPulmonaryVeins);
Aaron Bray's avatar
Aaron Bray committed
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
}

//--------------------------------------------------------------------------------------------------
/// \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()
{
410
  /// \todo figure out how to save an initial healthy/baseline patient definition and another current patient definition
411
  //Experimentally determined  
412
  double tidalVolumeBaseline_L = GetTidalVolume(VolumeUnit::L);
413
  //Calculated
414
415
416
  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;
Aaron Bray's avatar
Aaron Bray committed
417
418
419
420
421
422

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

  std::stringstream ss;
423
  ss << typeString << "Patient tidal volume = " << tidalVolumeBaseline_L << " L.";
Aaron Bray's avatar
Aaron Bray committed
424
  Info(ss);
425
  ss << typeString << "Patient inspiratory reserve volume = " << inspiratoryReserveVolume_L << " L.";
Aaron Bray's avatar
Aaron Bray committed
426
427
428
  Info(ss);

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

    //Update healthy patient volumes
434
435
436
437
    //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);
Aaron Bray's avatar
Aaron Bray committed
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
  }
}

//--------------------------------------------------------------------------------------------------
/// \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()
{
454
455
456
  CalculateWork();
  CalculateFatigue();

457
  UpdateChestWallCompliances();
Aaron Bray's avatar
Aaron Bray committed
458
459
460
461
462
463
464
  UpdateVolumes();
  UpdateResistances();
  UpdateAlveolarCompliances();
  UpdateInspiratoryExpiratoryRatio();
  UpdateDiffusion();
  UpdatePulmonaryCapillary();

Aaron Bray's avatar
Aaron Bray committed
465
466
467
468
469
  ProcessAerosolSubstances();
  Pneumothorax();  
  ConsciousRespiration();

  MechanicalVentilation();
470
  SupplementalOxygen();
Aaron Bray's avatar
Aaron Bray committed
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491

  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();
492
  
Aaron Bray's avatar
Aaron Bray committed
493
  // Calc the circuits
494
  m_Calculator->Process(RespirationCircuit, m_dt_s);
495
  
Aaron Bray's avatar
Aaron Bray committed
496
497
498
  //ModifyPleuralVolume();
  SEGasCompartmentGraph& RespirationGraph = m_data.GetCompartments().GetActiveRespiratoryGraph();
  SELiquidCompartmentGraph& AerosolGraph = m_data.GetCompartments().GetActiveAerosolGraph();
499
  
Aaron Bray's avatar
Aaron Bray committed
500
  // Transport substances
501
  m_GasTransporter->Transport(RespirationGraph, m_dt_s);
Aaron Bray's avatar
Aaron Bray committed
502
  if(m_AerosolMouth->HasSubstanceQuantities())
503
    m_AerosolTransporter->Transport(AerosolGraph, m_dt_s);
504
  
Aaron Bray's avatar
Aaron Bray committed
505
506
  //Update system data
  CalculateVitalSigns();
507

508
#ifdef DEBUG
509
  Debugging(RespirationCircuit);
510
#endif  
Aaron Bray's avatar
Aaron Bray committed
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
}

//--------------------------------------------------------------------------------------------------
/// \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();  
526
  m_Calculator->PostProcess(RespirationCircuit);  
Aaron Bray's avatar
Aaron Bray committed
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
}

//--------------------------------------------------------------------------------------------------
/// \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
614
615
    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();
Aaron Bray's avatar
Aaron Bray committed
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
    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
640
641
    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();
Aaron Bray's avatar
Aaron Bray committed
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
    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)
668
    if (subQ->GetSubstance().GetState() == eSubstance_State::Liquid)
Aaron Bray's avatar
Aaron Bray committed
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
    {
      // 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);  

696
  // Change resistances due to deposition
697
  double dTracheaResistance = m_MouthToCarina->GetNextResistance(PressureTimePerVolumeUnit::cmH2O_s_Per_L)*carinaResistanceModifier;
698
699
  double dLeftAlveoliResistance = m_LeftAnatomicDeadSpaceToLeftAlveolarDeadSpace->GetNextResistance(PressureTimePerVolumeUnit::cmH2O_s_Per_L)*leftAlveoliResistanceModifier;
  double dRightAlveoliResistance = m_RightAnatomicDeadSpaceToRightAlveolarDeadSpace->GetNextResistance(PressureTimePerVolumeUnit::cmH2O_s_Per_L)*rightAlveoliResistanceModifier;
Aaron Bray's avatar
Aaron Bray committed
700

701
702
  double dLeftBronchiResistance = m_CarinaToLeftAnatomicDeadSpace->GetNextResistance(PressureTimePerVolumeUnit::cmH2O_s_Per_L);
  double dRightBronchiResistance = m_CarinaToRightAnatomicDeadSpace->GetNextResistance(PressureTimePerVolumeUnit::cmH2O_s_Per_L);
Aaron Bray's avatar
Aaron Bray committed
703
704
  dLeftBronchiResistance *= leftDeadSpaceResistanceModifier * combinedLeftBronchodilationEffects;
  dRightBronchiResistance *= rightDeadSpaceResistanceModifier * combinedRightBronchodilationEffects;
705
706
  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);
Aaron Bray's avatar
Aaron Bray committed
707

708
709
710
  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);
711
712
  m_LeftAnatomicDeadSpaceToLeftAlveolarDeadSpace->GetNextResistance().SetValue(dLeftAlveoliResistance, PressureTimePerVolumeUnit::cmH2O_s_Per_L);
  m_RightAnatomicDeadSpaceToRightAlveolarDeadSpace->GetNextResistance().SetValue(dRightAlveoliResistance, PressureTimePerVolumeUnit::cmH2O_s_Per_L);
Aaron Bray's avatar
Aaron Bray committed
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
}

//--------------------------------------------------------------------------------------------------
/// \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
729
    m_data.SetAirwayMode(eAirwayMode::MechanicalVentilator);
730
731
732
733
734
735
736
737
738
739

    //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())
Aaron Bray's avatar
Aaron Bray committed
740
    {
741
742
743
744
      for (auto s : m_Environment->GetSubstanceQuantities())
      {
        m_MechanicalVentilatorConnection->GetSubstanceQuantity(s->GetSubstance())->GetVolumeFraction().Set(s->GetVolumeFraction());
      }
Aaron Bray's avatar
Aaron Bray committed
745
    }
746
    else
Aaron Bray's avatar
Aaron Bray committed
747
    {
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
      //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();
Aaron Bray's avatar
Aaron Bray committed
765

766
767
768
    //Reset the substance quantities at the connection
    for (SELiquidSubstanceQuantity* subQ : m_MechanicalVentilatorAerosolConnection->GetSubstanceQuantities())
      subQ->SetToZero();
Aaron Bray's avatar
Aaron Bray committed
769

770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
    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);
      }
Aaron Bray's avatar
Aaron Bray committed
785
786
    }

787
788
    //Apply the instantaneous flow ********************************************
    if (mv->HasFlow())
Aaron Bray's avatar
Aaron Bray committed
789
    {
790
791
792
793
794
795
      m_ConnectionToMouth->GetNextFlowSource().Set(mv->GetFlow());
      //It may or may not be there
      if (!m_ConnectionToMouth->HasFlowSource())
      {
        m_data.GetCircuits().GetRespiratoryAndMechanicalVentilatorCircuit().StateChange();
      }
Aaron Bray's avatar
Aaron Bray committed
796
    }
797
    else
Aaron Bray's avatar
Aaron Bray committed
798
    {
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
      //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);
Aaron Bray's avatar
Aaron Bray committed
817
818
    }
  }
819
  else if (m_data.GetAirwayMode() == eAirwayMode::MechanicalVentilator)
Aaron Bray's avatar
Aaron Bray committed
820
821
  {
    // Was just turned off
822
    m_data.SetAirwayMode(eAirwayMode::Free);
823
824
825
826
827
    if (m_ConnectionToMouth->HasNextFlowSource())
    {
      m_ConnectionToMouth->GetNextFlowSource().Invalidate();
      m_data.GetCircuits().GetRespiratoryAndMechanicalVentilatorCircuit().StateChange();
    }
Aaron Bray's avatar
Aaron Bray committed
828
829
830
  }
}

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

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

847
848
849
850
851
852
853
  // 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())
854
  {
855
    case eSupplementalOxygen_Device::None:
856
    {
857
858
859
      m_data.SetAirwayMode(eAirwayMode::Free);
      m_data.GetActions().GetPatientActions().RemoveSupplementalOxygen();
      return;
860
    }
861
    case eSupplementalOxygen_Device::NasalCannula:
862
    {
863
864
865
866
867
868
869
870
871
872
873
      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;
874
    }
875
    case eSupplementalOxygen_Device::NonRebreatherMask:
876
    {
877
878
879
880
881
882
883
884
885
886
887
      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;
888
    }
889
    case eSupplementalOxygen_Device::SimpleMask:
890
    {
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
      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;
907
908
909
910
911
912
913
    }
  }

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

  //Decrement volume from the tank
919
920
921
922
923
924
925
926
927
928
  //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.
Aaron Bray's avatar
Aaron Bray committed
929
      m_data.GetEvents().SetEvent(eEvent::SupplementalOxygenBottleExhausted, true, m_data.GetSimulationTime());
930
    }
Aaron Bray's avatar
Aaron Bray committed
931
932
    else
      m_data.GetEvents().SetEvent(eEvent::SupplementalOxygenBottleExhausted, false, m_data.GetSimulationTime());
933
  }
Aaron Bray's avatar
Aaron Bray committed
934
935
936
  else
    m_data.GetEvents().SetEvent(eEvent::SupplementalOxygenBottleExhausted, false, m_data.GetSimulationTime());

937
938
  //Nonrebreather mask works differently with the bag and doesn't have a pressure source for the tank
  if (so->GetDevice() != eSupplementalOxygen_Device::NonRebreatherMask)
939
940
941
942
943
944
945
946
  {
    //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
947
    OxygenInlet->GetNextResistance().SetValue(resistance_cmH2O_s_Per_L, PressureTimePerVolumeUnit::cmH2O_s_Per_L);
948
949
950
951
  }
  else
  {
    //Handle nonrebreather mask bag
952
953
954
    SEFluidCircuitNode* NonRebreatherMaskBag = RespirationCircuit->GetNode(pulse::NonRebreatherMaskNode::NonRebreatherMaskBag);
    SEFluidCircuitPath* NonRebreatherMaskReservoirValve = RespirationCircuit->GetPath(pulse::NonRebreatherMaskPath::NonRebreatherMaskReservoirValve);

955
    double flowOut_L_Per_min = 0.0;
956
    if (NonRebreatherMaskReservoirValve->HasNextFlow())
957
    {
958
      flowOut_L_Per_min = NonRebreatherMaskReservoirValve->GetNextFlow(VolumePerTimeUnit::L_Per_min);
959
    }
960

961
962
    double bagVolume_L = NonRebreatherMaskBag->GetNextVolume(VolumeUnit::L);
    bagVolume_L = bagVolume_L + (flow_L_Per_min - flowOut_L_Per_min) * m_dt_min;
963
964
965
    if (bagVolume_L < 0.0)
    {
      bagVolume_L = 0.0;
966

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

      /// \event Supplemental Oxygen: The nonrebreather mask is empty. Oxygen may need to be provided at a faster rate.
Aaron Bray's avatar
Aaron Bray committed
971
      m_data.GetEvents().SetEvent(eEvent::NonRebreatherMaskOxygenBagEmpty, true, m_data.GetSimulationTime());
972
    }
Aaron Bray's avatar
Aaron Bray committed
973
    else
974
    {
Aaron Bray's avatar
Aaron Bray committed
975
976
977
978
979
      if (bagVolume_L > 1.0)
      {
        bagVolume_L = 1.0;
      }
      m_data.GetEvents().SetEvent(eEvent::NonRebreatherMaskOxygenBagEmpty, false, m_data.GetSimulationTime());
980
981
    }
    
982
    NonRebreatherMaskBag->GetNextVolume().SetValue(bagVolume_L, VolumeUnit::L);
983
  }
984
985
}

Aaron Bray's avatar
Aaron Bray committed
986
987
988
989
990
991
992
993
994
995
996
997
998
//--------------------------------------------------------------------------------------------------
/// \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()
{
Jeff Webb's avatar
Jeff Webb committed
999
  /// \event Patient: Start of exhale/inhale
1000
1001
1002
1003
  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());
Jeff Webb's avatar
Jeff Webb committed
1004

Aaron Bray's avatar
Aaron Bray committed
1005
1006
1007
  m_BreathingCycleTime_s += m_dt_s;

  //Keep a running average of the Arterial Partial Pressures  
1008
1009
  m_ArterialO2RunningAverage_mmHg->Sample(m_AortaO2->GetPartialPressure(PressureUnit::mmHg));
  m_ArterialCO2RunningAverage_mmHg->Sample(m_AortaCO2->GetPartialPressure(PressureUnit::mmHg));
Aaron Bray's avatar
Aaron Bray committed
1010
  //Reset at start of cardiac cycle 
1011
  if (m_data.GetEvents().IsEventActive(eEvent::StartOfCardiacCycle))
Aaron Bray's avatar
Aaron Bray committed
1012
  {
1013
1014
    m_ArterialO2PartialPressure_mmHg = m_ArterialO2RunningAverage_mmHg->Value();
    m_ArterialCO2PartialPressure_mmHg = m_ArterialCO2RunningAverage_mmHg->Value();
Aaron Bray's avatar
Aaron Bray committed
1015

1016
1017
    m_ArterialO2RunningAverage_mmHg->Clear();
    m_ArterialCO2RunningAverage_mmHg->Clear();
Aaron Bray's avatar
Aaron Bray committed
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
  }

#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 
1039
      //pressure (attempting to reach a specific volume) at the end of the user specified period
Aaron Bray's avatar
Aaron Bray committed
1040
1041
1042
1043
1044
1045
1046
      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
1047
      m_VentilationFrequency_Per_min = m_data.GetCurrentPatient().GetRespirationRateBaseline(FrequencyUnit::Per_min);
1048
      m_BreathingCycleTime_s = 60.0 / m_VentilationFrequency_Per_min + 2.0 * m_dt_s;
Aaron Bray's avatar
Aaron Bray committed
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
    }
  }
  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
1069
      if (m_data.GetEvents().IsEventActive(eEvent::CardiacArrest))
Aaron Bray's avatar
Aaron Bray committed
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
      {
        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);
1080
      double BMR_W = m_data.GetCurrentPatient().GetBasalMetabolicRate(PowerUnit::W);
Aaron Bray's avatar
Aaron Bray committed
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
      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;
      }
1102
      else if (Drugs.GetSedationLevel().GetValue() > 0.15 && std::abs(m_data.GetCurrentPatient().GetRespirationRateBaseline(FrequencyUnit::Per_min) + DrugRRChange_Per_min) < 1.0)
Aaron Bray's avatar
Aaron Bray committed
1103
1104
1105
1106
1107
1108
1109
      {
        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
1110
1111
      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
Aaron Bray's avatar
Aaron Bray committed
1112

1113
1114
1115
      //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    
Aaron Bray's avatar
Aaron Bray committed
1116
      double tunedVolumeMetabolicSlope = 0.2; //Tuned fractional increase of the tidal volume due to increased metabolic rate
1117
      metabolicModifier = 1.0 + tunedVolumeMetabolicSlope * (metabolicFraction - 1.0);
Aaron Bray's avatar
Aaron Bray committed
1118
1119
1120
      dTargetAlveolarVentilation_L_Per_min *= metabolicModifier;

      //Only move it part way to where it wants to be.
1121
      //If it stays there, it will get there, just more controlled/slowly.
Aaron Bray's avatar
Aaron Bray committed
1122
1123
      //This is needed so we don't overshoot and introduce low frequency oscillations into the system (i.e. overdamped).
      double targetCO2PP_mmHg = 40.0;
1124
      double changeFraction = 1.0;
Aaron Bray's avatar
Aaron Bray committed
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
      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);
1136
1137
1138
1139
1140

#ifdef TUNING
      changeFraction = 1.0;
#endif

Aaron Bray's avatar
Aaron Bray committed
1141
1142
1143
1144
1145
      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
1146
      double dTargetPulmonaryVentilation_L_Per_min = dTargetAlveolarVentilation_L_Per_min;
Aaron Bray's avatar
Aaron Bray committed
1147
1148
1149
1150
1151
1152

      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;
1153
        m_data.GetEvents().SetEvent(eEvent::MaximumPulmonaryVentilationRate, true, m_data.GetSimulationTime());
Aaron Bray's avatar
Aaron Bray committed
1154
1155
      }

1156
      if (dTargetPulmonaryVentilation_L_Per_min < dMaximumPulmonaryVentilationRate && m_data.GetEvents().IsEventActive(eEvent::MaximumPulmonaryVentilationRate))
Aaron Bray's avatar
Aaron Bray committed
1157
      {
1158
        m_data.GetEvents().SetEvent(eEvent::MaximumPulmonaryVentilationRate, false, m_data.GetSimulationTime());
Aaron Bray's avatar
Aaron Bray committed
1159
1160
1161
1162
1163
1164
      }

      //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
1165
      dTargetTidalVolume_L *= cardiacArrestEffect * NMBModifier;
Aaron Bray's avatar
Aaron Bray committed
1166
1167
1168
1169
1170
1171

      //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
1172
      double dHalfVitalCapacity_L = m_data.GetCurrentPatient().GetVitalCapacity(VolumeUnit::L) / 2;
Aaron Bray's avatar
Aaron Bray committed
1173
1174
1175
      dTargetTidalVolume_L = MIN(dTargetTidalVolume_L, dHalfVitalCapacity_L);

      //Map the Target Tidal Volume to the Driver
1176
      double TargetVolume_L = m_data.GetCurrentPatient().GetFunctionalResidualCapacity(VolumeUnit::L) + dTargetTidalVolume_L;
Aaron Bray's avatar
Aaron Bray committed
1177
1178
1179
      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);
1180

Aaron Bray's avatar
Aaron Bray committed
1181
1182
1183
1184
1185
1186

      //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;
1187
        m_NotBreathing = true;
Aaron Bray's avatar
Aaron Bray committed
1188
1189
1190
1191
1192
1193
      }
      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;
1194
        m_NotBreathing = false;
Aaron Bray's avatar
Aaron Bray committed
1195
1196
      }

1197
      m_VentilationFrequency_Per_min = LIMIT(m_VentilationFrequency_Per_min, 0.0, dMaximumPulmonaryVentilationRate / dHalfVitalCapacity_L);
Aaron Bray's avatar
Aaron Bray committed
1198
1199
1200
1201
1202
1203

      //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)
      {
1204
        double dRespirationRateBaseline_Per_min = m_data.GetCurrentPatient().GetRespirationRateBaseline(FrequencyUnit::Per_min);
Aaron Bray's avatar
Aaron Bray committed
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
        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
1215
      m_VentilationToTidalVolumeSlope = 40.0;
Aaron Bray's avatar
Aaron Bray committed
1216
1217
1218
      m_VentilationFrequency_Per_min = 16.0;
#endif

1219
1220
      SetBreathCycleFractions();

Aaron Bray's avatar
Aaron Bray committed
1221
1222
1223
1224
1225
      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;
1226

1227
      m_data.GetEvents().SetEvent(eEvent::StartOfInhale, true, m_data.GetSimulationTime()); ///\todo rename "StartOfConsiousInhale"
Aaron Bray's avatar
Aaron Bray committed
1228
1229
1230
    }

    //Run the driver based on the waveform -------------------------------------------------------------------------------
1231
1232
1233
    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
Aaron Bray's avatar
Aaron Bray committed
1234

1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
    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