Last active
March 23, 2022 23:27
-
-
Save icook/899fcc7187d2d0ae5f57 to your computer and use it in GitHub Desktop.
A little class (and python notebook) to calculate and plot the concentration and buildup of Vyvanse in your body over time.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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)) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment