Module pk2Comp

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

# In[51]:


# Import commands
from scipy.stats import gamma
import numpy as np
import matplotlib.pyplot as plt
#%matplotlib inline
from scipy.integrate import odeint 
import math as math
import os
import pandas as pd

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

    def __init__ (self, numParam = 4, Flow = 1/60, Vp = 0.05, Visf = 0.15, PS = 1/60):
        
        """Initializes the model with default parameter values for flow, Vp, Visf, and PS.
        Parameters
        ----------      
        numParam: int
            numParam is the number of parameters you want to optimize for the model. Defaults to 4.
            
        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.      
        """
        
        # Declare Variables for initial conditions
        self.numParam = numParam
        self.flow = Flow
        self.Vp = Vp
        self.Visf = Visf
        self.PS = PS
        self.sol = []
        self.Cp0 = 0 # Initial concentration of tracer in plasma
        self.Cisf0 = 0 # Initial concentration of tracer in interstitial space
        self.tmax = 10 #Time in seconds
        self.dt = 0.1 #Time step
        self.a = 2. # Alpha for gamma distribution
        self.rv = gamma(self.a, loc = 2, scale = 0.55) #input function
        self.sol = []
        self.Mass_plasma = [] #mass of tracer in plasma
        self.Mass_isf = []
        self.Q = []
        
        # Define the time array
        self.time = np.arange(0, self.tmax + self.dt, self.dt)
        
    
    # Derivative function
    def derivs(self, curr_vals, time):
        """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
        -------
        dCp_dt : double[]
            contains the derivative of concentration in plasma with respect to time.
        dCisf_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
        Cp, Cisf = curr_vals

        # Define value of input function Cin
        Cin = self.rv.pdf(time)    

        # Right-hand side of odes, which are used to computer the derivative
        dCp_dt = (self.flow/self.Vp)*(Cin - Cp) + (self.PS/self.Vp)*(Cisf - Cp)
        dCisf_dt = (self.PS/self.Visf)*(Cp - Cisf)

        return dCp_dt, dCisf_dt

    def main(self):
        """Main function to solve ODEs"""
        # Store the initial values in a list
        init = [self.Cp0, self.Cisf0]

        # Solve the odes with odeint
        self.sol = odeint(self.derivs, init, self.time)

        self.Mass_plasma = self.Vp * self.sol[:,0] #mass of tracer in plasma
        self.Mass_isf = self.Visf * self.sol[:,1] #mass of tracer in isf
        #Tp = Vp/(flow + PS) # mean transit time
        #E = 1 - np.exp(-PS/flow) #extraction fraction
        self.Q = self.Mass_plasma + self.Mass_isf

        #print('The mean transit time is ' + str(Tp))
        #print('The extraction fraction is ' + str(E))

    def getPlot(self):
        """Plots the solution of the solved ODEs.
        
        Attributes
        ----------      
        sol : double[]
            contains the solutions of our ODE functions. 
        """

        # Plot the results using the values stored in the solution variable, "sol"
        # Plot Cp using the "0" element from the solution
        plt.figure(1)
        plt.plot(self.time, self.rv.pdf(self.time), color = 'blue', label = 'Input Function')
        plt.plot(self.time, self.sol[:,0],color="green", label = 'Cp')

        # Plot Cisf using the "1" element from the solution
        plt.plot(self.time, self.sol[:,1],color="purple", label = 'Cisf')
        plt.xlabel('Time [s]')
        plt.ylabel('Concentration [mM]')
        plt.legend(loc = 'best')
        plt.grid()

        # Plot mass of tracer using the "2" element from the solution
        plt.figure(2)
        plt.plot(self.time, self.Mass_plasma,color="red", label = 'Plasma')
        
        # Plot mass of tracer in tissue using the "3" element from the solution
        plt.plot(self.time, self.Mass_isf,color="black", label = 'Interstitial Space')
        plt.plot(self.time, self.Q, color="blue", label = 'Total mass')
        plt.xlabel('Time [s]')
        plt.ylabel('Mass [mg]')
        plt.legend(loc = 'best')
        plt.grid()

        print('Cp at 10 sec is ' + str(self.sol[100,0]))
        print('Cisf at 10 sec is ' + str(self.sol[100,1]))


# In[52]:


test = pk2Comp(4,.2,.2,.7,1/60)
test.main()
test.getPlot()


# In[ ]:

Classes

class pk2Comp (numParam=4, Flow=0.016666666666666666, Vp=0.05, Visf=0.15, PS=0.016666666666666666)

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


numParam: int numParam is the number of parameters you want to optimize for the model. Defaults to 4.

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.

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

    def __init__ (self, numParam = 4, Flow = 1/60, Vp = 0.05, Visf = 0.15, PS = 1/60):
        
        """Initializes the model with default parameter values for flow, Vp, Visf, and PS.
        Parameters
        ----------      
        numParam: int
            numParam is the number of parameters you want to optimize for the model. Defaults to 4.
            
        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.      
        """
        
        # Declare Variables for initial conditions
        self.numParam = numParam
        self.flow = Flow
        self.Vp = Vp
        self.Visf = Visf
        self.PS = PS
        self.sol = []
        self.Cp0 = 0 # Initial concentration of tracer in plasma
        self.Cisf0 = 0 # Initial concentration of tracer in interstitial space
        self.tmax = 10 #Time in seconds
        self.dt = 0.1 #Time step
        self.a = 2. # Alpha for gamma distribution
        self.rv = gamma(self.a, loc = 2, scale = 0.55) #input function
        self.sol = []
        self.Mass_plasma = [] #mass of tracer in plasma
        self.Mass_isf = []
        self.Q = []
        
        # Define the time array
        self.time = np.arange(0, self.tmax + self.dt, self.dt)
        
    
    # Derivative function
    def derivs(self, curr_vals, time):
        """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
        -------
        dCp_dt : double[]
            contains the derivative of concentration in plasma with respect to time.
        dCisf_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
        Cp, Cisf = curr_vals

        # Define value of input function Cin
        Cin = self.rv.pdf(time)    

        # Right-hand side of odes, which are used to computer the derivative
        dCp_dt = (self.flow/self.Vp)*(Cin - Cp) + (self.PS/self.Vp)*(Cisf - Cp)
        dCisf_dt = (self.PS/self.Visf)*(Cp - Cisf)

        return dCp_dt, dCisf_dt

    def main(self):
        """Main function to solve ODEs"""
        # Store the initial values in a list
        init = [self.Cp0, self.Cisf0]

        # Solve the odes with odeint
        self.sol = odeint(self.derivs, init, self.time)

        self.Mass_plasma = self.Vp * self.sol[:,0] #mass of tracer in plasma
        self.Mass_isf = self.Visf * self.sol[:,1] #mass of tracer in isf
        #Tp = Vp/(flow + PS) # mean transit time
        #E = 1 - np.exp(-PS/flow) #extraction fraction
        self.Q = self.Mass_plasma + self.Mass_isf

        #print('The mean transit time is ' + str(Tp))
        #print('The extraction fraction is ' + str(E))

    def getPlot(self):
        """Plots the solution of the solved ODEs.
        
        Attributes
        ----------      
        sol : double[]
            contains the solutions of our ODE functions. 
        """

        # Plot the results using the values stored in the solution variable, "sol"
        # Plot Cp using the "0" element from the solution
        plt.figure(1)
        plt.plot(self.time, self.rv.pdf(self.time), color = 'blue', label = 'Input Function')
        plt.plot(self.time, self.sol[:,0],color="green", label = 'Cp')

        # Plot Cisf using the "1" element from the solution
        plt.plot(self.time, self.sol[:,1],color="purple", label = 'Cisf')
        plt.xlabel('Time [s]')
        plt.ylabel('Concentration [mM]')
        plt.legend(loc = 'best')
        plt.grid()

        # Plot mass of tracer using the "2" element from the solution
        plt.figure(2)
        plt.plot(self.time, self.Mass_plasma,color="red", label = 'Plasma')
        
        # Plot mass of tracer in tissue using the "3" element from the solution
        plt.plot(self.time, self.Mass_isf,color="black", label = 'Interstitial Space')
        plt.plot(self.time, self.Q, color="blue", label = 'Total mass')
        plt.xlabel('Time [s]')
        plt.ylabel('Mass [mg]')
        plt.legend(loc = 'best')
        plt.grid()

        print('Cp at 10 sec is ' + str(self.sol[100,0]))
        print('Cisf at 10 sec is ' + str(self.sol[100,1]))

Methods

def derivs(self, curr_vals, time)

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

dCp_dt : double[]
contains the derivative of concentration in plasma with respect to time.
dCisf_dt : double[]
contains the derivative of concentration in interstitial fluid with respect to time.
Expand source code
def derivs(self, curr_vals, time):
    """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
    -------
    dCp_dt : double[]
        contains the derivative of concentration in plasma with respect to time.
    dCisf_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
    Cp, Cisf = curr_vals

    # Define value of input function Cin
    Cin = self.rv.pdf(time)    

    # Right-hand side of odes, which are used to computer the derivative
    dCp_dt = (self.flow/self.Vp)*(Cin - Cp) + (self.PS/self.Vp)*(Cisf - Cp)
    dCisf_dt = (self.PS/self.Visf)*(Cp - Cisf)

    return dCp_dt, dCisf_dt
def getPlot(self)

Plots the solution of the solved ODEs.

Attributes

sol : double[] contains the solutions of our ODE functions.

Expand source code
def getPlot(self):
    """Plots the solution of the solved ODEs.
    
    Attributes
    ----------      
    sol : double[]
        contains the solutions of our ODE functions. 
    """

    # Plot the results using the values stored in the solution variable, "sol"
    # Plot Cp using the "0" element from the solution
    plt.figure(1)
    plt.plot(self.time, self.rv.pdf(self.time), color = 'blue', label = 'Input Function')
    plt.plot(self.time, self.sol[:,0],color="green", label = 'Cp')

    # Plot Cisf using the "1" element from the solution
    plt.plot(self.time, self.sol[:,1],color="purple", label = 'Cisf')
    plt.xlabel('Time [s]')
    plt.ylabel('Concentration [mM]')
    plt.legend(loc = 'best')
    plt.grid()

    # Plot mass of tracer using the "2" element from the solution
    plt.figure(2)
    plt.plot(self.time, self.Mass_plasma,color="red", label = 'Plasma')
    
    # Plot mass of tracer in tissue using the "3" element from the solution
    plt.plot(self.time, self.Mass_isf,color="black", label = 'Interstitial Space')
    plt.plot(self.time, self.Q, color="blue", label = 'Total mass')
    plt.xlabel('Time [s]')
    plt.ylabel('Mass [mg]')
    plt.legend(loc = 'best')
    plt.grid()

    print('Cp at 10 sec is ' + str(self.sol[100,0]))
    print('Cisf at 10 sec is ' + str(self.sol[100,1]))
def main(self)

Main function to solve ODEs

Expand source code
def main(self):
    """Main function to solve ODEs"""
    # Store the initial values in a list
    init = [self.Cp0, self.Cisf0]

    # Solve the odes with odeint
    self.sol = odeint(self.derivs, init, self.time)

    self.Mass_plasma = self.Vp * self.sol[:,0] #mass of tracer in plasma
    self.Mass_isf = self.Visf * self.sol[:,1] #mass of tracer in isf
    #Tp = Vp/(flow + PS) # mean transit time
    #E = 1 - np.exp(-PS/flow) #extraction fraction
    self.Q = self.Mass_plasma + self.Mass_isf