diff --git a/pk_optimizer/pkOptimizer.py b/pk_optimizer/pkOptimizer.py
deleted file mode 100644
index d5423c7929acc838efb9376d36aa0fe079e73fd6..0000000000000000000000000000000000000000
--- a/pk_optimizer/pkOptimizer.py
+++ /dev/null
@@ -1,216 +0,0 @@
-#!/usr/bin/env python
-# coding: utf-8
-
-# In[1]:
-
-
-from scipy.stats import gamma
-from scipy.integrate import odeint 
-from scipy.optimize import fmin as fmin
-
-import os
-import csv
-import re
-import math as math
-import numpy as np
-import matplotlib.pyplot as plt
-
-class pkOptimizer:
-    """The pkOptimizer object is an optimizer for parameters in pk models."""
-    
-    def __init__ (self, wd):#, Flow = 1/60, Vp = 0.05, Visf = 0.15, PS = 1/60):
-        """Initializes the model with initial guess parameter values for flow, Vp, Visf, and PS.
-        Parameters
-        ----------      
-        wd : string
-            wd is the working directory path
-        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.
-        time : double[]
-            list of all timepoints
-        aorta : double[]
-            concentration of tracer in aorta (input function)
-        myo : double[]
-            concentration of tracer in myocardial tissue (Cisf)
-        opt : float[]
-            opt is the optimized parameters to fit curve
-        """
-        self.wd = wd
-        #self.Flow = Flow
-        #self.Vp = Vp
-        #self.Visf = Visf
-        #self.PS = PS
-        self.time = []
-        self.aorta = []
-        self.myo = []
-        self.opt = []
-        
-    def getData(self, filename):
-        """Imports data from all .csv files in directory.
-        Parameters
-        ----------  
-        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 (Cisf)
-        
-        Returns
-        -------
-        time : double[]
-            list of all timepoints
-        aorta : double[]
-            concentration of tracer in aorta (input function)
-        myo : double[]
-            concentration of tracer in myocardial tissue (Cisf)
-        """ 
-        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('\d+[.]+\d+|\d+').findall(data[i+1][0])[0]))
-            self.aorta.append(float(re.compile('\d+[.]+\d+|\d+').findall(data[i+1][1])[0]))
-            self.myo.append(float(re.compile('\d+[.]+\d+|\d+').findall(data[i+1][2])[0]))
-            
-        return self.time, self.aorta, self.myo
-
-    #gamma_var distribution curve
-    def gamma_var(self, t = np.arange(0,25), ymax = 10, tmax = 10, alpha = 2, delay = 5):
-        """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 = t
-        ymax = ymax
-        tmax = tmax
-        alpha = alpha
-        delay = delay
-
-        y = np.zeros(np.size(t)); #preallocate output
-
-        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
-
-    #gamma_var_error
-    def MSE(self, param = [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(param) < 1:
-            ymax = 10; tmax = 10; alpha = 2,; delay = 5
-        elif len(param) < 2:
-            ymax = param[0]; tmax = 10; alpha = 2,; delay = 5
-        elif len(param) < 3:
-            ymax = param[0]; tmax = param[1]; alpha = 2,; delay = 5
-        elif len(param) < 4:
-            ymax = param[0]; tmax = param[1]; alpha = param[2],; delay = 5
-        else:
-            # Mean squared error (MSE) between data and gamma variate with given parameters
-            ymax = param[0]
-            tmax = param[1]
-            alpha = param[2]
-            delay = param[3]
-        
-        MSE = 0
-
-        if tmax <= 0 or ymax <= 10 or delay < 0 or alpha<0 or alpha>1000 or tmax>1000:
-            MSE = 1000000; #just return a big number
-            return MSE
-
-        else:
-            model_vals = self.gamma_var(self.time, ymax, tmax, alpha, delay);
-
-            for i in range(len(self.myo)):
-                MSE = (self.aorta[i]-model_vals[i])**2 + MSE
-            MSE = MSE / len(self.myo)           
-            return round(MSE,3)
-
-    def inputFuncFit(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
-        self.opt = fmin(self.MSE, initGuesses, maxiter = 1000)          
-        return self.opt.round(2)
-    
-    def plotOpt(self):
-        """Plots the original data to the fitted curve."""
-        plt.plot(self.time, self.aorta, 'bo', label='data')
-        #plt.plot(t, y, 'b-', label='data')
-        
-        y = self.gamma_var(np.arange(0,25,0.01), self.opt[0], self.opt[1], self.opt[2], self.opt[3])
-        
-        plt.plot(np.arange(0,25,0.01), y, 'r-')
-
-