Module pk_one_comp

The pk1Comp object is a one compartment PK model that outputs graphs of mass of tracer over time.

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

# In[1]:
"""The pk1Comp object is a one compartment PK
    model that outputs graphs of mass of tracer over time."""

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_one_comp:

    """The pk1Comp object is a one compartment PK
    model that outputs graphs of mass 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 :
            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

        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)
        """
        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)
            return y

        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, init_guesses):
        """Uses fmin algorithm (Nelder-Mead simplex algorithm) to
        minimize loss function (MSE) of input function.
        Parameters
        ----------
        init_guesses : 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, init_guesses, 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
        -------
        dc_dt : double[]
            contains the derivative of concentrations with respect to time.
        """

        # Define value of input function conc_in
        conc_in = self.gamma_var(time, self.ymax, self.tmax, self.alpha, self.delay)
        # Unpack the current values of the variables we wish to "update" from the curr_vals list
        myo_vals = curr_vals

        # Right-hand side of odes, which are used to computer the derivative
        dc_dt = self.flow*(conc_in - myo_vals)/self.vol_plasma
        return dc_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.vol_plasma = guess[1]
        self.baseline = guess[2]

        mse = 0

        if self.flow <= 0 or self.flow >= 25 or self.vol_plasma > 100 or self.vol_plasma < 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], 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] - sol.y[0][i])**2 + mse

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

    def output_func_fit(self, init_guesses):
        """Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize
        loss function (MSE) of output function.
        Parameters
        ----------
        init_guesses : 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, init_guesses, maxiter=10000)

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

        return opt1.round(4)

    def main(self):
        """Run this for demonstration of what pk_one_comp can do"""
        # 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])
        print(opt)
        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([.03102, .462, self.myo[0]])
        print(opt2)
        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_one_comp (wd=WindowsPath('C:/Users/Ethan/OneDrive - Michigan State University/MSU/Classwork/Computational Modeling/Models/Data'), filename='CTPERF005_stress.csv')

The pk1Comp object is a one compartment PK model that outputs graphs of mass 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 : 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/(mLmin). 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/(gmin). 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_one_comp:

    """The pk1Comp object is a one compartment PK
    model that outputs graphs of mass 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 :
            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

        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)
        """
        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)
            return y

        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, init_guesses):
        """Uses fmin algorithm (Nelder-Mead simplex algorithm) to
        minimize loss function (MSE) of input function.
        Parameters
        ----------
        init_guesses : 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, init_guesses, 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
        -------
        dc_dt : double[]
            contains the derivative of concentrations with respect to time.
        """

        # Define value of input function conc_in
        conc_in = self.gamma_var(time, self.ymax, self.tmax, self.alpha, self.delay)
        # Unpack the current values of the variables we wish to "update" from the curr_vals list
        myo_vals = curr_vals

        # Right-hand side of odes, which are used to computer the derivative
        dc_dt = self.flow*(conc_in - myo_vals)/self.vol_plasma
        return dc_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.vol_plasma = guess[1]
        self.baseline = guess[2]

        mse = 0

        if self.flow <= 0 or self.flow >= 25 or self.vol_plasma > 100 or self.vol_plasma < 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], 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] - sol.y[0][i])**2 + mse

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

    def output_func_fit(self, init_guesses):
        """Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize
        loss function (MSE) of output function.
        Parameters
        ----------
        init_guesses : 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, init_guesses, maxiter=10000)

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

        return opt1.round(4)

    def main(self):
        """Run this for demonstration of what pk_one_comp can do"""
        # 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])
        print(opt)
        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([.03102, .462, self.myo[0]])
        print(opt2)
        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

dc_dt : double[]
contains the derivative of concentrations 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
    -------
    dc_dt : double[]
        contains the derivative of concentrations with respect to time.
    """

    # Define value of input function conc_in
    conc_in = self.gamma_var(time, self.ymax, self.tmax, self.alpha, self.delay)
    # Unpack the current values of the variables we wish to "update" from the curr_vals list
    myo_vals = curr_vals

    # Right-hand side of odes, which are used to computer the derivative
    dc_dt = self.flow*(conc_in - myo_vals)/self.vol_plasma
    return dc_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)
        return y

    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

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

    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)
    """
    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, init_guesses)

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


init_guesses : 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, init_guesses):
    """Uses fmin algorithm (Nelder-Mead simplex algorithm) to
    minimize loss function (MSE) of input function.
    Parameters
    ----------
    init_guesses : 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, init_guesses, 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)

Run this for demonstration of what pk_one_comp can do

Expand source code
def main(self):
    """Run this for demonstration of what pk_one_comp can do"""
    # 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])
    print(opt)
    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([.03102, .462, self.myo[0]])
    print(opt2)
    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, init_guesses)

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


init_guesses : 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, init_guesses):
    """Uses fmin algorithm (Nelder-Mead simplex algorithm) to minimize
    loss function (MSE) of output function.
    Parameters
    ----------
    init_guesses : 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, init_guesses, maxiter=10000)

    self.flow = opt1[0].round(4)
    self.vol_plasma = 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.vol_plasma = guess[1]
    self.baseline = guess[2]

    mse = 0

    if self.flow <= 0 or self.flow >= 25 or self.vol_plasma > 100 or self.vol_plasma < 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], 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] - sol.y[0][i])**2 + mse

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