Module pk_two_comp

The pk2Comp object is a two compartment PK model that outputs graphs of concentration of tracer over time.

Expand source code
#!/usr/bin/env python
# coding: utf-8

# In[1]:


"""The pk2Comp object is a two compartment PK model
    that outputs graphs of concentration of tracer over time.
"""
#!/usr/bin/env python
# coding: utf-8

# In[1]:
import pathlib
import os
import csv
import re
import math
import numpy as np
import matplotlib.pyplot as plt

from scipy.integrate import solve_ivp
from scipy.optimize import fmin

class pk_two_comp:
    """The pk2Comp object is a two compartment PK model
    that outputs graphs of concentration of tracer over time.
    """

    def __init__(self, wd=pathlib.Path('Data').absolute(), filename='CTPERF005_stress.csv'):
        """Initializes the model with default parameter values for flow, Vp, Visf, and PS.
        Parameters
        ----------
        wd : path
            Absolute path name to current directory. Defaults to ./Data
        filename : String
            Name of the data file you'd like to access

        Attributes
        -----------
        time : double[]
            list of all timepoints
        aorta : double[]
            concentration of tracer in aorta (input function)
        myo : double[]
            concentration of tracer in myocardial tissue (conc_isf)
        Flow : double
            Flow is the flow of plasma through the blood vessel in mL/(mL*min). Defaults to 1/60.
        Vp : double
            Vp is the volume of plasma in mL. Defaults to 0.05.
        Visf : double
            Visf is the volume of interstitial fluid in mL. Defaults to 0.15.
        PS : double
            PS is the permeability-surface area constant in mL/(g*min). Defaults to 1/60.
        ymax : int
            Magnitude of Gamma-var peak.
        tmax : double
            Time of which highest peak in Gamma-var appears
        alpha : double
            Scale factor
        delay : double
            Delay to start Gamma-var curve.

        """
        # Subject Data
        if os.path.basename(os.path.normpath(pathlib.Path().absolute())) != 'Data':
            self.wd = pathlib.Path('Data').absolute()
        else:
            self.wd = wd

        if not isinstance(filename, str):
            raise ValueError('Filename must be a string')

        self.filename = filename
        self.time = []
        self.aorta = []
        self.myo = []

        # Declare Variables for initial conditions
        # Compartment variables to be fitted
        self.flow = 1/60
        self.visf = 0.15
        self.baseline = 60

        # Other Compartmental Modelvariables
        self.perm_surf = 0.35
        self.vol_plasma = 0.10

        # Solved ode
        self.sol = []

        # Gamma variables
        self.ymax = 250
        self.tmax = 6.5
        self.alpha = 2.5
        self.delay = 0

        self.deriv_sol = np.array([])
        self.fit_myo = np.array([])

    def get_data(self, filename):
        """Imports data from all .csv files in directory.
        Parameters
        ----------
        filename : string
            Name of the data file you'd like to access
        wd : str
            wd is the working directory path

        Attributes
        ----------
        time : double[]
            list of all timepoints
        aorta : double[]
            concentration of tracer in aorta (input function)
        myo : double[]
            concentration of tracer in myocardial tissue (conc_isf)

        Returns
        -------
        time : double[]
            list of all timepoints
        aorta : double[]
            concentration of tracer in aorta (input function)
        myo : double[]
            concentration of tracer in myocardial tissue (conc_isf)
        """
        self.time = []
        self.aorta = []
        self.myo = []

        os.chdir(self.wd)
        # File not found error
        if not os.path.isfile(filename):
            raise ValueError(
                "Input file does not exist: {0}. I'll quit now.".format(filename))

        data = list(csv.reader(open(filename), delimiter='\t'))

        for i in range(12):
            self.time.append(
                float(re.compile(r'\d+[.]+\d+|\d+').findall(data[i+1][0])[0]))
            self.aorta.append(
                float(re.compile(r'\d+[.]+\d+|\d+').findall(data[i+1][1])[0]))
            self.myo.append(
                float(re.compile(r'\d+[.]+\d+|\d+').findall(data[i+1][2])[0]))

        return self.time, self.aorta, self.myo

    # gamma_var distribution curve
    def gamma_var(self, time=np.arange(0, 25), ymax=10, tmax=10, alpha=2, delay=0):
        """Creates a gamma variate probability density function with given alpha,
        location, and scale values.
        Parameters
        ----------
        t : double[]
            array of timepoints
        ymax : double
            peak y value of gamma distribution
        tmax : double
            location of 50th percentile of function
        alpha : double
            scale parameter
        delay : double
            time delay of which to start gamma distribution

        Returns
        -------
        y : double[]
            probability density function of your gamma variate.
        """
        # Following Madsen 1992 simplified parameterization for gamma variate
        t = time
        self.ymax = ymax
        self.tmax = tmax
        self.alpha = alpha
        self.delay = delay

        y = np.zeros(np.size(t))  # preallocate output

        # For odeint, checks if t input is array or float
        if isinstance(t, (list, np.ndarray)):
            for i in range(np.size(y)):
                if t[i] < delay:
                    y[i] = 0
                else:
                    y[i] = round((ymax*tmax**(-alpha)*math.exp(alpha))*(t[i]-delay)
                                 ** alpha*math.exp(-alpha*(t[i]-delay)/tmax), 3)
        else:
            y = round((ymax*tmax**(-alpha)*math.exp(alpha))*(t-delay)
                      ** alpha*math.exp(-alpha*(t-delay)/tmax), 3)
        return y

    # gamma_var_error
    def input_mse(self, guess=[10, 10, 2, 5]):
        """Calculates Mean squared error (MSE) between data and
        gamma variate with given parameters values.
        Parameters
        ----------
        param : ndarray[]
            time : double[]
                array of timepoints
            ymax : double
                peak y value of gamma distribution
            tmax : double
                location of 50th percentile of function
            alpha : double
                scale parameter
            delay : double
                time delay of which to start gamma distribution

        Returns
        -------
        MSE : double
            Mean squared error
        """
        if len(guess) < 1:
            self.ymax = 10
            self.tmax = 10
            self.alpha = 2
            self.delay = 5
        elif len(guess) < 2:
            self.ymax = guess[0]
            self.tmax = 10
            self.alpha = 2
            self.delay = 5
        elif len(guess) < 3:
            self.ymax = guess[0]
            self.tmax = guess[1]
            self.alpha = 2
            self.delay = 5
        elif len(guess) < 4:
            self.ymax = guess[0]
            self.tmax = guess[1]
            self.alpha = guess[2]
            self.delay = 5
        else:
            # Mean squared error (MSE) between data and gamma variate with given parameters
            self.ymax = guess[0]
            self.tmax = guess[1]
            self.alpha = guess[2]
            self.delay = guess[3]

        mse = 0

        if self.tmax <= 0 or self.ymax <= 10 or self.delay < 0 or self.alpha < 0 \
            or self.alpha > 1000 or self.tmax > 1000:
            mse = 1000000  # just return a big number

        else:
            model_vals = self.gamma_var(
                self.time, self.ymax, self.tmax, self.alpha, self.delay)

            for i in range(len(self.aorta)):
                mse = (self.aorta[i] - model_vals[i])**2 + mse
            mse = mse / len(self.aorta)
        return round(mse, 3)

    def input_func_fit(self, initGuesses):
        """Uses fmin algorithm (Nelder-Mead simplex algorithm) to
        minimize loss function (MSE) of input function.
        Parameters
        ----------
        initGuesses : ndarray[]
            Array of initial guesses containing:
                time : double[]
                    array of timepoints
                ymax : double
                    peak y value of gamma distribution
                tmax : double
                    location of 50th percentile of function
                alpha : double
                    scale parameter
                delay : double
                    time delay of which to start gamma distribution
        Returns
        -------
        opt : double[]
            optimized parameters
        """
        # Mean squared error (MSE) between data and gamma variate with given parameters
        opt = fmin(self.input_mse, initGuesses, maxiter=1000)

        self.ymax = opt[0]
        self.tmax = opt[1]
        self.alpha = opt[2]
        self.delay = opt[3]

        return opt.round(2)

    # Derivative function
    def derivs(self, time, curr_vals):
        """Finds derivatives of ODEs.

        Parameters
        ----------
        curr_vals : double[]
            curr_vals it he current values of the variables we wish to
            "update" from the curr_vals list.

        time : double[]
            time is our time array from 0 to tmax with timestep dt.

        Returns
        -------
        dconc_plasma_dt : double[]
            contains the derivative of concentration in plasma with respect to time.
        dconc_isf_dt : double[]
            contains the derivative of concentration in interstitial fluid with respect to time.
        """
        # Unpack the current values of the variables we wish to "update" from the curr_vals list
        conc_plasma, conc_isf = curr_vals

        # Define value of input function conc_in
        conc_in = self.gamma_var(time, self.ymax, self.tmax,
                                 self.alpha, self.delay)

        # Right-hand side of odes, which are used to computer the derivative
        dconc_plasma_dt = (self.flow/self.vol_plasma)*(conc_in - conc_plasma) + \
            (self.perm_surf/self.vol_plasma)*(conc_isf - conc_plasma)
        dconc_isf_dt = (self.perm_surf/self.visf)*(conc_plasma - conc_isf)
        return dconc_plasma_dt, dconc_isf_dt

    def output_mse(self, guess):
        """Calculates Mean squared error (MSE) between data and
        output derivs with given parameters values.
        Parameters
        ----------
        guess : ndarray[]

            Flow : double
                Flow is the flow of plasma through the blood vessel in mL/(mL*min).
                Defaults to 1/60.

            Vp : double
                Vp is the volume of plasma in mL. Defaults to 0.05.

            Visf : double
                Visf is the volume of interstitial fluid in mL. Defaults to 0.15.

            PS : double
                PS is the permeability-surface area constant in mL/(g*min). Defaults to 1/60.
        Returns
        -------
        MSE : double
            Mean squared error
        """
        self.flow = guess[0]
        self.visf = guess[1]
        self.baseline = guess[2]

        mse = 0

        if self.flow <= 0 or self.flow >= 25 or self.visf > 100 \
            or self.visf < 0 or self.baseline > 150 or self.baseline < 0:
            mse = 100000  # just return a big number

        else:
            sol = solve_ivp(self.derivs, [0, 30], [0, 0], t_eval=self.time)
            MBF = sol.y[0] + sol.y[1]

            temp = np.asarray(self.myo) - self.baseline

            for i in range(len(self.myo)):
                mse = (temp[i] - MBF[i])**2 + mse

            mse = mse / len(self.myo)
        return mse

    def output_func_fit(self, initGuesses):
        """Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize
        loss function (MSE) of output function.
        Parameters
        ----------
        initGuesses : ndarray[]
            Array of initial guesses containing:
                flow : double
                    Flow is the flow of plasma through the blood vessel in mL/(mL*min).
                    Defaults to 1/60.

            Vp : double
                Vp is the volume of plasma in mL. Defaults to 0.05.

            Visf : double
                Visf is the volume of interstitial fluid in mL. Defaults to 0.15.

            PS : double
                PS is the permeability-surface area constant in mL/(g*min). Defaults to 1/60.
        Returns
        -------
        opt : double[]
            optimized parameters
        """
        # Mean squared error (MSE) between data and gamma variate with given parameters
        opt1 = fmin(self.output_mse, initGuesses, maxiter=10000)

        self.flow = opt1[0].round(4)
        self.visf = opt1[1].round(4)
        self.baseline = opt1[2].round(4)

        return opt1.round(4)

    def main(self):

        # Gets data from file
        self.get_data(self.filename)

        # Plots original data
        plt.plot(self.time, self.aorta, 'bo')
        plt.plot(self.time, self.myo, 'ro')

        # Fit gamma_var input function and plots it
        opt = self.input_func_fit([250, 7, 4, 0])
        plt.plot(np.arange(0, 25, 0.01), self.gamma_var(np.arange(0, 25, 0.01),
                                                        opt[0], opt[1], opt[2], opt[3]), 'k-')

        # Fit uptake function and plot it
        opt2 = self.output_func_fit([.011418, .62, self.myo[0]])
        self.deriv_sol = solve_ivp(self.derivs, [0, 30],
                                   [0, 0], t_eval=self.time)
        self.fit_myo = self.deriv_sol.y[0] + self.deriv_sol.y[1]
        plt.plot(self.time, self.fit_myo + self.baseline, 'm-')
        

Classes

class pk_two_comp (wd=WindowsPath('C:/Users/Ethan/OneDrive - Michigan State University/MSU/Classwork/Computational Modeling/Models/Data'), filename='CTPERF005_stress.csv')

The pk2Comp object is a two compartment PK model that outputs graphs of concentration of tracer over time.

Initializes the model with default parameter values for flow, Vp, Visf, and PS. Parameters


wd : path
Absolute path name to current directory. Defaults to ./Data
filename : String
Name of the data file you'd like to access

Attributes

time : double[]
list of all timepoints
aorta : double[]
concentration of tracer in aorta (input function)
myo : double[]
concentration of tracer in myocardial tissue (conc_isf)
Flow : double
Flow is the flow of plasma through the blood vessel in mL/(mL*min). Defaults to 1/60.
Vp : double
Vp is the volume of plasma in mL. Defaults to 0.05.
Visf : double
Visf is the volume of interstitial fluid in mL. Defaults to 0.15.
PS : double
PS is the permeability-surface area constant in mL/(g*min). Defaults to 1/60.
ymax : int
Magnitude of Gamma-var peak.
tmax : double
Time of which highest peak in Gamma-var appears
alpha : double
Scale factor
delay : double
Delay to start Gamma-var curve.
Expand source code
class pk_two_comp:
    """The pk2Comp object is a two compartment PK model
    that outputs graphs of concentration of tracer over time.
    """

    def __init__(self, wd=pathlib.Path('Data').absolute(), filename='CTPERF005_stress.csv'):
        """Initializes the model with default parameter values for flow, Vp, Visf, and PS.
        Parameters
        ----------
        wd : path
            Absolute path name to current directory. Defaults to ./Data
        filename : String
            Name of the data file you'd like to access

        Attributes
        -----------
        time : double[]
            list of all timepoints
        aorta : double[]
            concentration of tracer in aorta (input function)
        myo : double[]
            concentration of tracer in myocardial tissue (conc_isf)
        Flow : double
            Flow is the flow of plasma through the blood vessel in mL/(mL*min). Defaults to 1/60.
        Vp : double
            Vp is the volume of plasma in mL. Defaults to 0.05.
        Visf : double
            Visf is the volume of interstitial fluid in mL. Defaults to 0.15.
        PS : double
            PS is the permeability-surface area constant in mL/(g*min). Defaults to 1/60.
        ymax : int
            Magnitude of Gamma-var peak.
        tmax : double
            Time of which highest peak in Gamma-var appears
        alpha : double
            Scale factor
        delay : double
            Delay to start Gamma-var curve.

        """
        # Subject Data
        if os.path.basename(os.path.normpath(pathlib.Path().absolute())) != 'Data':
            self.wd = pathlib.Path('Data').absolute()
        else:
            self.wd = wd

        if not isinstance(filename, str):
            raise ValueError('Filename must be a string')

        self.filename = filename
        self.time = []
        self.aorta = []
        self.myo = []

        # Declare Variables for initial conditions
        # Compartment variables to be fitted
        self.flow = 1/60
        self.visf = 0.15
        self.baseline = 60

        # Other Compartmental Modelvariables
        self.perm_surf = 0.35
        self.vol_plasma = 0.10

        # Solved ode
        self.sol = []

        # Gamma variables
        self.ymax = 250
        self.tmax = 6.5
        self.alpha = 2.5
        self.delay = 0

        self.deriv_sol = np.array([])
        self.fit_myo = np.array([])

    def get_data(self, filename):
        """Imports data from all .csv files in directory.
        Parameters
        ----------
        filename : string
            Name of the data file you'd like to access
        wd : str
            wd is the working directory path

        Attributes
        ----------
        time : double[]
            list of all timepoints
        aorta : double[]
            concentration of tracer in aorta (input function)
        myo : double[]
            concentration of tracer in myocardial tissue (conc_isf)

        Returns
        -------
        time : double[]
            list of all timepoints
        aorta : double[]
            concentration of tracer in aorta (input function)
        myo : double[]
            concentration of tracer in myocardial tissue (conc_isf)
        """
        self.time = []
        self.aorta = []
        self.myo = []

        os.chdir(self.wd)
        # File not found error
        if not os.path.isfile(filename):
            raise ValueError(
                "Input file does not exist: {0}. I'll quit now.".format(filename))

        data = list(csv.reader(open(filename), delimiter='\t'))

        for i in range(12):
            self.time.append(
                float(re.compile(r'\d+[.]+\d+|\d+').findall(data[i+1][0])[0]))
            self.aorta.append(
                float(re.compile(r'\d+[.]+\d+|\d+').findall(data[i+1][1])[0]))
            self.myo.append(
                float(re.compile(r'\d+[.]+\d+|\d+').findall(data[i+1][2])[0]))

        return self.time, self.aorta, self.myo

    # gamma_var distribution curve
    def gamma_var(self, time=np.arange(0, 25), ymax=10, tmax=10, alpha=2, delay=0):
        """Creates a gamma variate probability density function with given alpha,
        location, and scale values.
        Parameters
        ----------
        t : double[]
            array of timepoints
        ymax : double
            peak y value of gamma distribution
        tmax : double
            location of 50th percentile of function
        alpha : double
            scale parameter
        delay : double
            time delay of which to start gamma distribution

        Returns
        -------
        y : double[]
            probability density function of your gamma variate.
        """
        # Following Madsen 1992 simplified parameterization for gamma variate
        t = time
        self.ymax = ymax
        self.tmax = tmax
        self.alpha = alpha
        self.delay = delay

        y = np.zeros(np.size(t))  # preallocate output

        # For odeint, checks if t input is array or float
        if isinstance(t, (list, np.ndarray)):
            for i in range(np.size(y)):
                if t[i] < delay:
                    y[i] = 0
                else:
                    y[i] = round((ymax*tmax**(-alpha)*math.exp(alpha))*(t[i]-delay)
                                 ** alpha*math.exp(-alpha*(t[i]-delay)/tmax), 3)
        else:
            y = round((ymax*tmax**(-alpha)*math.exp(alpha))*(t-delay)
                      ** alpha*math.exp(-alpha*(t-delay)/tmax), 3)
        return y

    # gamma_var_error
    def input_mse(self, guess=[10, 10, 2, 5]):
        """Calculates Mean squared error (MSE) between data and
        gamma variate with given parameters values.
        Parameters
        ----------
        param : ndarray[]
            time : double[]
                array of timepoints
            ymax : double
                peak y value of gamma distribution
            tmax : double
                location of 50th percentile of function
            alpha : double
                scale parameter
            delay : double
                time delay of which to start gamma distribution

        Returns
        -------
        MSE : double
            Mean squared error
        """
        if len(guess) < 1:
            self.ymax = 10
            self.tmax = 10
            self.alpha = 2
            self.delay = 5
        elif len(guess) < 2:
            self.ymax = guess[0]
            self.tmax = 10
            self.alpha = 2
            self.delay = 5
        elif len(guess) < 3:
            self.ymax = guess[0]
            self.tmax = guess[1]
            self.alpha = 2
            self.delay = 5
        elif len(guess) < 4:
            self.ymax = guess[0]
            self.tmax = guess[1]
            self.alpha = guess[2]
            self.delay = 5
        else:
            # Mean squared error (MSE) between data and gamma variate with given parameters
            self.ymax = guess[0]
            self.tmax = guess[1]
            self.alpha = guess[2]
            self.delay = guess[3]

        mse = 0

        if self.tmax <= 0 or self.ymax <= 10 or self.delay < 0 or self.alpha < 0 \
            or self.alpha > 1000 or self.tmax > 1000:
            mse = 1000000  # just return a big number

        else:
            model_vals = self.gamma_var(
                self.time, self.ymax, self.tmax, self.alpha, self.delay)

            for i in range(len(self.aorta)):
                mse = (self.aorta[i] - model_vals[i])**2 + mse
            mse = mse / len(self.aorta)
        return round(mse, 3)

    def input_func_fit(self, initGuesses):
        """Uses fmin algorithm (Nelder-Mead simplex algorithm) to
        minimize loss function (MSE) of input function.
        Parameters
        ----------
        initGuesses : ndarray[]
            Array of initial guesses containing:
                time : double[]
                    array of timepoints
                ymax : double
                    peak y value of gamma distribution
                tmax : double
                    location of 50th percentile of function
                alpha : double
                    scale parameter
                delay : double
                    time delay of which to start gamma distribution
        Returns
        -------
        opt : double[]
            optimized parameters
        """
        # Mean squared error (MSE) between data and gamma variate with given parameters
        opt = fmin(self.input_mse, initGuesses, maxiter=1000)

        self.ymax = opt[0]
        self.tmax = opt[1]
        self.alpha = opt[2]
        self.delay = opt[3]

        return opt.round(2)

    # Derivative function
    def derivs(self, time, curr_vals):
        """Finds derivatives of ODEs.

        Parameters
        ----------
        curr_vals : double[]
            curr_vals it he current values of the variables we wish to
            "update" from the curr_vals list.

        time : double[]
            time is our time array from 0 to tmax with timestep dt.

        Returns
        -------
        dconc_plasma_dt : double[]
            contains the derivative of concentration in plasma with respect to time.
        dconc_isf_dt : double[]
            contains the derivative of concentration in interstitial fluid with respect to time.
        """
        # Unpack the current values of the variables we wish to "update" from the curr_vals list
        conc_plasma, conc_isf = curr_vals

        # Define value of input function conc_in
        conc_in = self.gamma_var(time, self.ymax, self.tmax,
                                 self.alpha, self.delay)

        # Right-hand side of odes, which are used to computer the derivative
        dconc_plasma_dt = (self.flow/self.vol_plasma)*(conc_in - conc_plasma) + \
            (self.perm_surf/self.vol_plasma)*(conc_isf - conc_plasma)
        dconc_isf_dt = (self.perm_surf/self.visf)*(conc_plasma - conc_isf)
        return dconc_plasma_dt, dconc_isf_dt

    def output_mse(self, guess):
        """Calculates Mean squared error (MSE) between data and
        output derivs with given parameters values.
        Parameters
        ----------
        guess : ndarray[]

            Flow : double
                Flow is the flow of plasma through the blood vessel in mL/(mL*min).
                Defaults to 1/60.

            Vp : double
                Vp is the volume of plasma in mL. Defaults to 0.05.

            Visf : double
                Visf is the volume of interstitial fluid in mL. Defaults to 0.15.

            PS : double
                PS is the permeability-surface area constant in mL/(g*min). Defaults to 1/60.
        Returns
        -------
        MSE : double
            Mean squared error
        """
        self.flow = guess[0]
        self.visf = guess[1]
        self.baseline = guess[2]

        mse = 0

        if self.flow <= 0 or self.flow >= 25 or self.visf > 100 \
            or self.visf < 0 or self.baseline > 150 or self.baseline < 0:
            mse = 100000  # just return a big number

        else:
            sol = solve_ivp(self.derivs, [0, 30], [0, 0], t_eval=self.time)
            MBF = sol.y[0] + sol.y[1]

            temp = np.asarray(self.myo) - self.baseline

            for i in range(len(self.myo)):
                mse = (temp[i] - MBF[i])**2 + mse

            mse = mse / len(self.myo)
        return mse

    def output_func_fit(self, initGuesses):
        """Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize
        loss function (MSE) of output function.
        Parameters
        ----------
        initGuesses : ndarray[]
            Array of initial guesses containing:
                flow : double
                    Flow is the flow of plasma through the blood vessel in mL/(mL*min).
                    Defaults to 1/60.

            Vp : double
                Vp is the volume of plasma in mL. Defaults to 0.05.

            Visf : double
                Visf is the volume of interstitial fluid in mL. Defaults to 0.15.

            PS : double
                PS is the permeability-surface area constant in mL/(g*min). Defaults to 1/60.
        Returns
        -------
        opt : double[]
            optimized parameters
        """
        # Mean squared error (MSE) between data and gamma variate with given parameters
        opt1 = fmin(self.output_mse, initGuesses, maxiter=10000)

        self.flow = opt1[0].round(4)
        self.visf = opt1[1].round(4)
        self.baseline = opt1[2].round(4)

        return opt1.round(4)

    def main(self):

        # Gets data from file
        self.get_data(self.filename)

        # Plots original data
        plt.plot(self.time, self.aorta, 'bo')
        plt.plot(self.time, self.myo, 'ro')

        # Fit gamma_var input function and plots it
        opt = self.input_func_fit([250, 7, 4, 0])
        plt.plot(np.arange(0, 25, 0.01), self.gamma_var(np.arange(0, 25, 0.01),
                                                        opt[0], opt[1], opt[2], opt[3]), 'k-')

        # Fit uptake function and plot it
        opt2 = self.output_func_fit([.011418, .62, self.myo[0]])
        self.deriv_sol = solve_ivp(self.derivs, [0, 30],
                                   [0, 0], t_eval=self.time)
        self.fit_myo = self.deriv_sol.y[0] + self.deriv_sol.y[1]
        plt.plot(self.time, self.fit_myo + self.baseline, 'm-')

Methods

def derivs(self, time, curr_vals)

Finds derivatives of ODEs.

Parameters

curr_vals : double[]
curr_vals it he current values of the variables we wish to "update" from the curr_vals list.
time : double[]
time is our time array from 0 to tmax with timestep dt.

Returns

dconc_plasma_dt : double[]
contains the derivative of concentration in plasma with respect to time.
dconc_isf_dt : double[]
contains the derivative of concentration in interstitial fluid with respect to time.
Expand source code
def derivs(self, time, curr_vals):
    """Finds derivatives of ODEs.

    Parameters
    ----------
    curr_vals : double[]
        curr_vals it he current values of the variables we wish to
        "update" from the curr_vals list.

    time : double[]
        time is our time array from 0 to tmax with timestep dt.

    Returns
    -------
    dconc_plasma_dt : double[]
        contains the derivative of concentration in plasma with respect to time.
    dconc_isf_dt : double[]
        contains the derivative of concentration in interstitial fluid with respect to time.
    """
    # Unpack the current values of the variables we wish to "update" from the curr_vals list
    conc_plasma, conc_isf = curr_vals

    # Define value of input function conc_in
    conc_in = self.gamma_var(time, self.ymax, self.tmax,
                             self.alpha, self.delay)

    # Right-hand side of odes, which are used to computer the derivative
    dconc_plasma_dt = (self.flow/self.vol_plasma)*(conc_in - conc_plasma) + \
        (self.perm_surf/self.vol_plasma)*(conc_isf - conc_plasma)
    dconc_isf_dt = (self.perm_surf/self.visf)*(conc_plasma - conc_isf)
    return dconc_plasma_dt, dconc_isf_dt
def gamma_var(self, time=array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]), ymax=10, tmax=10, alpha=2, delay=0)

Creates a gamma variate probability density function with given alpha, location, and scale values. Parameters


t : double[]
array of timepoints
ymax : double
peak y value of gamma distribution
tmax : double
location of 50th percentile of function
alpha : double
scale parameter
delay : double
time delay of which to start gamma distribution

Returns

y : double[]
probability density function of your gamma variate.
Expand source code
def gamma_var(self, time=np.arange(0, 25), ymax=10, tmax=10, alpha=2, delay=0):
    """Creates a gamma variate probability density function with given alpha,
    location, and scale values.
    Parameters
    ----------
    t : double[]
        array of timepoints
    ymax : double
        peak y value of gamma distribution
    tmax : double
        location of 50th percentile of function
    alpha : double
        scale parameter
    delay : double
        time delay of which to start gamma distribution

    Returns
    -------
    y : double[]
        probability density function of your gamma variate.
    """
    # Following Madsen 1992 simplified parameterization for gamma variate
    t = time
    self.ymax = ymax
    self.tmax = tmax
    self.alpha = alpha
    self.delay = delay

    y = np.zeros(np.size(t))  # preallocate output

    # For odeint, checks if t input is array or float
    if isinstance(t, (list, np.ndarray)):
        for i in range(np.size(y)):
            if t[i] < delay:
                y[i] = 0
            else:
                y[i] = round((ymax*tmax**(-alpha)*math.exp(alpha))*(t[i]-delay)
                             ** alpha*math.exp(-alpha*(t[i]-delay)/tmax), 3)
    else:
        y = round((ymax*tmax**(-alpha)*math.exp(alpha))*(t-delay)
                  ** alpha*math.exp(-alpha*(t-delay)/tmax), 3)
    return y
def get_data(self, filename)

Imports data from all .csv files in directory. Parameters


filename : string
Name of the data file you'd like to access
wd : str
wd is the working directory path

Attributes

time : double[]
list of all timepoints
aorta : double[]
concentration of tracer in aorta (input function)
myo : double[]
concentration of tracer in myocardial tissue (conc_isf)

Returns

time : double[]
list of all timepoints
aorta : double[]
concentration of tracer in aorta (input function)
myo : double[]
concentration of tracer in myocardial tissue (conc_isf)
Expand source code
def get_data(self, filename):
    """Imports data from all .csv files in directory.
    Parameters
    ----------
    filename : string
        Name of the data file you'd like to access
    wd : str
        wd is the working directory path

    Attributes
    ----------
    time : double[]
        list of all timepoints
    aorta : double[]
        concentration of tracer in aorta (input function)
    myo : double[]
        concentration of tracer in myocardial tissue (conc_isf)

    Returns
    -------
    time : double[]
        list of all timepoints
    aorta : double[]
        concentration of tracer in aorta (input function)
    myo : double[]
        concentration of tracer in myocardial tissue (conc_isf)
    """
    self.time = []
    self.aorta = []
    self.myo = []

    os.chdir(self.wd)
    # File not found error
    if not os.path.isfile(filename):
        raise ValueError(
            "Input file does not exist: {0}. I'll quit now.".format(filename))

    data = list(csv.reader(open(filename), delimiter='\t'))

    for i in range(12):
        self.time.append(
            float(re.compile(r'\d+[.]+\d+|\d+').findall(data[i+1][0])[0]))
        self.aorta.append(
            float(re.compile(r'\d+[.]+\d+|\d+').findall(data[i+1][1])[0]))
        self.myo.append(
            float(re.compile(r'\d+[.]+\d+|\d+').findall(data[i+1][2])[0]))

    return self.time, self.aorta, self.myo
def input_func_fit(self, initGuesses)

Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize loss function (MSE) of input function. Parameters


initGuesses : ndarray[]
Array of initial guesses containing: time : double[] array of timepoints ymax : double peak y value of gamma distribution tmax : double location of 50th percentile of function alpha : double scale parameter delay : double time delay of which to start gamma distribution

Returns

opt : double[]
optimized parameters
Expand source code
def input_func_fit(self, initGuesses):
    """Uses fmin algorithm (Nelder-Mead simplex algorithm) to
    minimize loss function (MSE) of input function.
    Parameters
    ----------
    initGuesses : ndarray[]
        Array of initial guesses containing:
            time : double[]
                array of timepoints
            ymax : double
                peak y value of gamma distribution
            tmax : double
                location of 50th percentile of function
            alpha : double
                scale parameter
            delay : double
                time delay of which to start gamma distribution
    Returns
    -------
    opt : double[]
        optimized parameters
    """
    # Mean squared error (MSE) between data and gamma variate with given parameters
    opt = fmin(self.input_mse, initGuesses, maxiter=1000)

    self.ymax = opt[0]
    self.tmax = opt[1]
    self.alpha = opt[2]
    self.delay = opt[3]

    return opt.round(2)
def input_mse(self, guess=[10, 10, 2, 5])

Calculates Mean squared error (MSE) between data and gamma variate with given parameters values. Parameters


param : ndarray[]
time : double[] array of timepoints ymax : double peak y value of gamma distribution tmax : double location of 50th percentile of function alpha : double scale parameter delay : double time delay of which to start gamma distribution

Returns

MSE : double
Mean squared error
Expand source code
def input_mse(self, guess=[10, 10, 2, 5]):
    """Calculates Mean squared error (MSE) between data and
    gamma variate with given parameters values.
    Parameters
    ----------
    param : ndarray[]
        time : double[]
            array of timepoints
        ymax : double
            peak y value of gamma distribution
        tmax : double
            location of 50th percentile of function
        alpha : double
            scale parameter
        delay : double
            time delay of which to start gamma distribution

    Returns
    -------
    MSE : double
        Mean squared error
    """
    if len(guess) < 1:
        self.ymax = 10
        self.tmax = 10
        self.alpha = 2
        self.delay = 5
    elif len(guess) < 2:
        self.ymax = guess[0]
        self.tmax = 10
        self.alpha = 2
        self.delay = 5
    elif len(guess) < 3:
        self.ymax = guess[0]
        self.tmax = guess[1]
        self.alpha = 2
        self.delay = 5
    elif len(guess) < 4:
        self.ymax = guess[0]
        self.tmax = guess[1]
        self.alpha = guess[2]
        self.delay = 5
    else:
        # Mean squared error (MSE) between data and gamma variate with given parameters
        self.ymax = guess[0]
        self.tmax = guess[1]
        self.alpha = guess[2]
        self.delay = guess[3]

    mse = 0

    if self.tmax <= 0 or self.ymax <= 10 or self.delay < 0 or self.alpha < 0 \
        or self.alpha > 1000 or self.tmax > 1000:
        mse = 1000000  # just return a big number

    else:
        model_vals = self.gamma_var(
            self.time, self.ymax, self.tmax, self.alpha, self.delay)

        for i in range(len(self.aorta)):
            mse = (self.aorta[i] - model_vals[i])**2 + mse
        mse = mse / len(self.aorta)
    return round(mse, 3)
def main(self)
Expand source code
def main(self):

    # Gets data from file
    self.get_data(self.filename)

    # Plots original data
    plt.plot(self.time, self.aorta, 'bo')
    plt.plot(self.time, self.myo, 'ro')

    # Fit gamma_var input function and plots it
    opt = self.input_func_fit([250, 7, 4, 0])
    plt.plot(np.arange(0, 25, 0.01), self.gamma_var(np.arange(0, 25, 0.01),
                                                    opt[0], opt[1], opt[2], opt[3]), 'k-')

    # Fit uptake function and plot it
    opt2 = self.output_func_fit([.011418, .62, self.myo[0]])
    self.deriv_sol = solve_ivp(self.derivs, [0, 30],
                               [0, 0], t_eval=self.time)
    self.fit_myo = self.deriv_sol.y[0] + self.deriv_sol.y[1]
    plt.plot(self.time, self.fit_myo + self.baseline, 'm-')
def output_func_fit(self, initGuesses)

Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize loss function (MSE) of output function. Parameters


initGuesses : ndarray[]

Array of initial guesses containing: flow : double Flow is the flow of plasma through the blood vessel in mL/(mL*min). Defaults to 1/60.

Vp : double Vp is the volume of plasma in mL. Defaults to 0.05.

Visf : double Visf is the volume of interstitial fluid in mL. Defaults to 0.15.

PS : double PS is the permeability-surface area constant in mL/(g*min). Defaults to 1/60.

Returns

opt : double[]
optimized parameters
Expand source code
def output_func_fit(self, initGuesses):
    """Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize
    loss function (MSE) of output function.
    Parameters
    ----------
    initGuesses : ndarray[]
        Array of initial guesses containing:
            flow : double
                Flow is the flow of plasma through the blood vessel in mL/(mL*min).
                Defaults to 1/60.

        Vp : double
            Vp is the volume of plasma in mL. Defaults to 0.05.

        Visf : double
            Visf is the volume of interstitial fluid in mL. Defaults to 0.15.

        PS : double
            PS is the permeability-surface area constant in mL/(g*min). Defaults to 1/60.
    Returns
    -------
    opt : double[]
        optimized parameters
    """
    # Mean squared error (MSE) between data and gamma variate with given parameters
    opt1 = fmin(self.output_mse, initGuesses, maxiter=10000)

    self.flow = opt1[0].round(4)
    self.visf = opt1[1].round(4)
    self.baseline = opt1[2].round(4)

    return opt1.round(4)
def output_mse(self, guess)

Calculates Mean squared error (MSE) between data and output derivs with given parameters values. Parameters


guess : ndarray[]

Flow : double Flow is the flow of plasma through the blood vessel in mL/(mL*min). Defaults to 1/60.

Vp : double Vp is the volume of plasma in mL. Defaults to 0.05.

Visf : double Visf is the volume of interstitial fluid in mL. Defaults to 0.15.

PS : double PS is the permeability-surface area constant in mL/(g*min). Defaults to 1/60.

Returns

MSE : double
Mean squared error
Expand source code
def output_mse(self, guess):
    """Calculates Mean squared error (MSE) between data and
    output derivs with given parameters values.
    Parameters
    ----------
    guess : ndarray[]

        Flow : double
            Flow is the flow of plasma through the blood vessel in mL/(mL*min).
            Defaults to 1/60.

        Vp : double
            Vp is the volume of plasma in mL. Defaults to 0.05.

        Visf : double
            Visf is the volume of interstitial fluid in mL. Defaults to 0.15.

        PS : double
            PS is the permeability-surface area constant in mL/(g*min). Defaults to 1/60.
    Returns
    -------
    MSE : double
        Mean squared error
    """
    self.flow = guess[0]
    self.visf = guess[1]
    self.baseline = guess[2]

    mse = 0

    if self.flow <= 0 or self.flow >= 25 or self.visf > 100 \
        or self.visf < 0 or self.baseline > 150 or self.baseline < 0:
        mse = 100000  # just return a big number

    else:
        sol = solve_ivp(self.derivs, [0, 30], [0, 0], t_eval=self.time)
        MBF = sol.y[0] + sol.y[1]

        temp = np.asarray(self.myo) - self.baseline

        for i in range(len(self.myo)):
            mse = (temp[i] - MBF[i])**2 + mse

        mse = mse / len(self.myo)
    return mse