Author: Daniele E. Schiavazzi
References:
Volume 33, Issue 3 March 2017 e02799
Date: June 28th, 2017
Objectives of this tutorial:
First, lets create a daData_multiple_Table tulip data object. This simply means that patient data will be stored by row where each patient occupies a separate column.
In [1]:
# Import tulip
import sys
sys.path.insert(0, '../tuliplib/tulipBin/py')
import numpy as np
import matplotlib.pyplot as plt
# Import UQ Library
import tulipUQ as uq
# Import Computational Model Library
import tulipCM as cm
# Import Data Library
import tulipDA as da
# Import Action Library
import tulipAC as ac
# READ DATASET
data = da.daData_multiple_Table()
data.readFromFile('heartFailure.dat')
print (np.loadtxt('heartFailure.dat',delimiter=',',dtype=str))
The columns in the above data table represent patients affected by diastolic left ventricular disfunction. Patients with normal, mild and severe left ventricular diastolic dysfunction are represented in columns 4, 3 and 2, repectively. As expected, the systemic indicators (Systemic Blood Pressures "SBP","DBP" and Cardiac Output "CO") are not significantly affected by this condition.
Let's start by creating a ODE circuit model. tulip provides a library of existing models, the odeNormalAdultSimplePA is an adult model with normal physiology and is characterized by a simple RC pulmonary compartment.
In [2]:
# CREATE ODE MODEL
ode = cm.odeNormalAdultSimplePA()
We now create a Runge-Kutta 4 ODE integrator and initialize it by providing the ODE model as weel as the integration time step and total number of heart cycles. Note that the ODE model knows about its own heart rate (one of the model parameters) and therefore the total number of time steps is automatically selected to consistently generated the same number of heart cycles.
In [3]:
# CREATE ODE INTEGRATOR
timeStep = 0.01
totalCycles = 10
rk4 = cm.odeIntegratorRK4(ode,timeStep,totalCycles)
An LPN (lumped parameter network) model is initialized by providing the selected ODE integrator.
In [4]:
# Create new LPN model
lpnModel = cm.cmLPNModel(rk4)
The cmLPNModel contains a data object that knows how to compute cost functions and likelihood. This object needs to be assigned before we can proceed. The second argument represents the specific column of the patient data file that we would like to use for optimization.
In [5]:
# ASSIGN DATA OBJECT TO MODEL
lpnModel.setData(data,0)
Once we are happy with the model and data, we need to decide what to do with them. One possibility is to do optimization. We first define some parameters and then initialize the acActionOPT_NM to perform optimization using the Nelder-Mead simplex algorithm.
The following parameters are used:
This algorithm has been modified from the implementation at this link.
References:
In [6]:
# SET OPTIMIZER PARAMETERS
# Total Number of iterations
totIterations = 3
# Convergence Tolerance
convTol = 1.0e-6
# Check Convergence every convUpdateIt iterations
convUpdateIt = 1
# Maximum Iterations
maxOptIt = 200
# Coefficient for Step increments
stepCoefficient = 0.1
# INIT ACTION
nm = ac.acActionOPT_NM(convTol,convUpdateIt,maxOptIt,stepCoefficient)
At this point, we assign the LPN model to the acAction object.
In [7]:
# ASSIGN MODEL TO ACTION
nm.setModel(lpnModel)
With this assignment we have completely defined the object hierarchy in tulip. Specifically, the model has a data member that, based on the model result, evaluates the likelihood or cost function and the optimizer knows about the model to run. This hierarchy is summarized in the picture below:
We set the initial guess to the default parameter set. This is a hardcoded set of parameters defined within each odeModel.
In [8]:
# SET INITIAL GUESS - DEFAULT PARAMETER VALUES
useStartingParameterFromFile = False
startFromCentre = False
startParameterFile = ''
nm.setInitialParamGuess(useStartingParameterFromFile,startFromCentre,startParameterFile)
We solve totIteration restarts using this loop. Note that the optimization is started by calling nm.go() and the iterations after the first will read the initial parameter guess by the file optParams.txt that contains the optimal set of parameters from the previous optimization.
In [9]:
# RUN RESTARTED NELDER-MEAD
for loopA in range(totIterations):
# PERFORM ACTION
nm.go()
# SET RESTART CONDITION
if(loopA == 0):
nm.setInitialPointFromFile(True)
nm.setInitialPointFile('optParams.txt')
print ('--- Iteration Completed!!')
Check the agreement between the patient data and the outputs of the optimized model. Let's open the outputTargets.out file that contains this information.
In [10]:
print (np.loadtxt('outputTargets.out',delimiter=',',dtype=str))
The allData.dat file has also been generated with the model outputs for all the time steps. Lets plot some pressures, flows and pressure-volume loops.
In [11]:
res = np.loadtxt('allData.dat')
ini = 500
# PLOT CONTENT
plt.figure(figsize=(20,5))
# Pressure Curves
plt.subplot(1,1,1)
plt.plot(res[ini:,11],res[ini:,19]/1333.3,'r-',lw=5,label='Left Ventricle') # Left Ventricular Pressure
plt.plot(res[ini:,11],res[ini:,17]/1333.3,'g-',lw=5,label='Left Atrium') # Left Atrial Pressure
plt.plot(res[ini:,11],res[ini:,18]/1333.3,'r--',lw=5,label='Right Ventricle') # Right Ventricular Pressure
plt.plot(res[ini:,11],res[ini:,16]/1333.3,'g--',lw=5,label='Right Atrium') # Right Atrial Pressure
plt.plot(res[ini:,11],res[ini:,8]/1333.3,'b-',lw=5,label='Aortic') # Aortic Pressure
plt.plot(res[ini:,11],res[ini:,5]/1333.3,'k-',lw=5,label='Pulmonary') # Pulmonary Pressure
plt.legend(fontsize=14)
plt.xlabel('Time [s]',fontsize=14)
plt.ylabel('Pressure [mmHg]',fontsize=14)
plt.tick_params(axis='both', which='major', labelsize=12)
plt.tick_params(axis='both', which='minor', labelsize=12)
plt.tight_layout()
plt.show()
In [12]:
# Flow Curves
plt.figure(figsize=(20,5))
plt.subplot(1,1,1)
plt.plot(res[ini:,11],res[ini:,4],'r-',lw=5,label='Rigth AV')
plt.plot(res[ini:,11],res[ini:,6],'g-',lw=5,label='Pulmonary')
plt.plot(res[ini:,11],res[ini:,7],'b-',lw=5,label='Left AV')
plt.plot(res[ini:,11],res[ini:,9],'k-',lw=5,label='Aortic')
plt.legend(fontsize=14)
plt.xlabel('Time [s]',fontsize=14)
plt.ylabel('Flow [cm$^3$/s]',fontsize=14)
plt.tick_params(axis='both', which='major', labelsize=12)
plt.tick_params(axis='both', which='minor', labelsize=12)
plt.tight_layout()
plt.show()
In [13]:
# Ventricular Pressure-Volume
plt.figure(figsize=(20,5))
plt.subplot(1,2,1)
plt.plot(res[ini:,3],res[ini:,19]/1333.3, lw=5) # Left Ventricular PV Loop
plt.xlabel('Left Ventricular Volume [cm$^3$]',fontsize=14)
plt.ylabel('Left Ventricular Pressure [mmHg]',fontsize=14)
plt.tick_params(axis='both', which='major', labelsize=12)
plt.tick_params(axis='both', which='minor', labelsize=12)
plt.subplot(1,2,2)
plt.plot(res[ini:,2],res[ini:,18]/1333.3, lw=5) # Right Ventricular PV Loop
plt.xlabel('Right Ventricular Volume [cm$^3$]',fontsize=14)
plt.ylabel('Right Ventricular Pressure [mmHg]',fontsize=14)
plt.tick_params(axis='both', which='major', labelsize=12)
plt.tick_params(axis='both', which='minor', labelsize=12)
plt.tight_layout()
plt.show()