A little class (and python notebook) to calculate and plot the concentration and buildup of Vyvanse in your body over time.

example.ipynb

simulation.py

import time
import datetime
import matplotlib.pyplot as plt
import numpy as np
import pandas
import pylab
from IPython.display import HTML
from scipy.signal import argrelextrema

pylab.rcParams['figure.figsize'] = 16, 8

PRODRUG_HALFLIFE = 3600
DEXTROAMPHETAMINE_HALFLIFE = 3600 * 11.5
print "Using {} seconds as halflife for metabolization of inactive component".format(PRODRUG_HALFLIFE)
print "Using {} seconds as halflife for dextroamphetamine".format(DEXTROAMPHETAMINE_HALFLIFE)


class Simulation(object):
    def __init__(self, **kwargs):
        self.__dict__.update(dict(dosage=70, interval=2, sample_interval=3600, duration=168, ingestion_interval=24, ingestion_start=8))
        self.__dict__.update(kwargs)
        self.interval = float(self.interval)
        now = datetime.datetime.now().replace(minute=0, hour=0, second=0, microsecond=0)
        self.start = now
        hours = 0
        prodrug_conc = 0.0
        dextroamphetamine_conc = 0.0
        ingestions = 0

        print "Configuration Details"
        print "Dosage: {}mg".format(self.dosage)
        print "Simulation Duration: {} days".format(self.duration/24)
        print "Taken every: {} hours".format(self.ingestion_interval)
        print "First taken at: {}\n".format(self.start.replace(hour=self.ingestion_start).strftime("%I:%M %p"))

        self.times = []
        self.unmet_conc = []
        self.met_conc = []

        self.ingestion_times = []
        self.ingestion_y = []
        for i in xrange(int(self.duration * 3600 / self.interval)):
            # At 8 am exactly
            ingestion = False
            if now == (self.start + datetime.timedelta(hours=self.ingestion_start + (self.ingestion_interval * ingestions))):
                ingestion = True
                ingestions += 1
                prodrug_conc += self.dosage
                #print "Took {}mg of Unmetabolized drug".format(dosage)

            # Reduce the active drug concentration by calculating half life decay
            dextroamphetamine_conc *= (0.5 ** (self.interval / DEXTROAMPHETAMINE_HALFLIFE))
            # Calculate the new concentration of prodrug in your system after the interval time
            new_prodrug_conc = prodrug_conc * (0.5 ** (self.interval / PRODRUG_HALFLIFE))
            # Increase the concentration of the active drug by the amount of prodrug lost
            dextroamphetamine_conc += prodrug_conc - new_prodrug_conc
            # Update the prodrug concentration
            prodrug_conc = new_prodrug_conc

            # If we took 70mg, flag it for the graph
            if ingestion:
                self.ingestion_times.append(hours)
                self.ingestion_y.append(prodrug_conc)

            # Print and record the current concentrations
            if time.mktime(now.timetuple()) % self.interval == 0:
                #print "At {} active drug concentration {:,.7f}. Unmetabolized concentration {:,.7f}.".format(now, dextroamphetamine_conc, prodrug_conc)
                self.times.append(hours)
                self.unmet_conc.append(prodrug_conc)
                self.met_conc.append(dextroamphetamine_conc)

            # increment the time
            now = now + datetime.timedelta(seconds=self.interval)
            hours = (now - self.start).total_seconds() / 3600

    def plot(self):
        ax = plt.subplot(111)
        ax.plot(self.times, self.unmet_conc, 'b', linewidth=2, label='Unmetabolized Concentration')
        ax.plot(self.times, self.met_conc, 'g', linewidth=2, label='Metabolized Concentration')
        ax.plot(self.ingestion_times, self.ingestion_y, 'ro', markersize=15, label='Ingestion Times')

        plt.title('Vyvanse dosage concentrations over {} days for {}mg taken '
                  'every {} hours starting at {}'
                  .format(int(self.duration / 24), self.dosage, self.ingestion_interval, self.ingestion_start))

        # Shrink the box to fit the legend
        box = ax.get_position()
        ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
        return (ax.legend(loc='center left', bbox_to_anchor=(1, 0.5)),
                )

    def concentration_at(self, target_hour):
        data = []
        for idx, hour in enumerate(self.times):
            if hour % 24 == target_hour:
                td = datetime.timedelta(hours=hour)
                date = self.start + td
                data.append(["{:,.1f}mg".format(self.met_conc[idx]),
                            "Day {}, {}".format(td.days, date.strftime("%I:%M %p"))])
        return HTML("<h3>Concentration at {}</h3>".format((self.start + datetime.timedelta(hours=target_hour)).strftime("%I:%M %p"))
                    + pandas.DataFrame.from_records(data, index=["Time"], columns=["Metabolized", "Time"]).to_html(escape=False))

    def min_max(self):
        data = []
        for idx in argrelextrema(np.array(self.met_conc), np.less)[0]:
            td = datetime.timedelta(hours=self.times[idx])
            date = self.start + td
            data.append(["Min", "{:,.1f}mg".format(self.met_conc[idx]),
                        "Day {}, {}".format(td.days, date.strftime("%I:%M %p")),
                        idx])

        for idx in argrelextrema(np.array(self.met_conc), np.greater)[0]:
            td = datetime.timedelta(hours=self.times[idx])
            date = self.start + td
            data.append(["Max", "{:,.1f}mg".format(self.met_conc[idx]),
                         "Day {}, {}".format(td.days, date.strftime("%I:%M %p")),
                         idx])
        data.sort(key=lambda s: s[3])
        data = [d[:3] for d in data]
        return HTML("<h3>Time of day for minimum and maximum drug concentrations</h3>"
                    + pandas.DataFrame.from_records(data, index=["Time"], columns=["Min/Max", "Metabolized", "Time"]).to_html(escape=False))