First try at reading Cp(t) then generating and plotting hysteresis curves.

Prettified from scratch_pad notebook, 2/17-18/2015.



In [35]:

    
import numpy as np

Definitions: time axis



In [36]:

    
time_step = 0.25
start_time = 0.0
duration = 480.0
n_rows = int(1+duration/time_step)
n_rows









    Out[36]:





1921

Definitions: PKPD parameters



In [37]:

    
ke0 = 5.91003E-03
nHill = 2.49464E+01
ec50 = 3.36280E+02

Open and read data file in CSV format



In [38]:

    
pwd









    Out[38]:





u'C:\\Users\\kevin\\Documents\\IPython_notebooks'

!xcopy C:\Users\kevin\Desktop\sync\hysteresis_phMRI\hysteresis_mapping_time_courses_20150217.csv .



In [39]:

    
with open('hysteresis_mapping_time_courses_20150217.csv', 'rb') as datafile:
    data = np.genfromtxt(datafile, delimiter=",", names=True, skip_header=3, usecols=(0,1,3,4,5,6,7,8,10,11,12,13,14,15))
Cp = data['Cp']
times = data['time_min']  # time, in minutes, from data file

WARNING This assumes data file time axis matches definitions above.

Alternatively, generate it from definitions above:

times = np.arange(start_time, duration+time_step, time_step)

Calculate Ce(t)

For now, try just our fake little Newton's method iteration. See this reference for alternative ODE solvers.

Note: I'm sure this step could be optimized, but it doesn't matter much at this stage. For one, given our previous efforts using Larry Bretthorst's Bayes Analysis Toolbox, scipy includes an F2Py tool to convert Fortran subroutines into callable Python functions.



In [44]:

    
Ce = np.zeros( n_rows )
for i in range(0,times.size-1) :
    Ce[i+1] = Ce[i] + (times[i+1]-times[i])*ke0*(Cp[i]-Ce[i])

Calculate effect(t)

Ignores scaling by Emax and shifting by e0, but operates on the whole array.



In [45]:

    
effect = np.zeros( n_rows )
effect = Ce**nHill/(ec50**nHill + Ce**nHill)

Plot Cp(t) and Ce(t), vs. EC50, and effect(t).

Alternatively:

fig, axes = plot.subplots(nrows=2, ncols=1, sharex=True, sharey=False)

See this page for instructions on plotting the effect on a separate, red, y-axis labeled on the right side of the graph.



In [46]:

    
import matplotlib.pyplot as plt
%matplotlib inline

figure = plt.figure(figsize=(5,6))  # figsize=(w,h) in inches
figure.clear

ax1 = figure.add_subplot(2,1,1)
ax1.plot(times,Cp,'black',linewidth=2)
ax1.plot(times,Ce,'black',linestyle='--',linewidth=2)
ax1.axhline(ec50,color='gray',linewidth=1)
ax1.axis((0,480,0,3200))
# ax1.set_xlabel('time (min)', fontsize=12)
ax1.set_ylabel('concentration (ng/ml)')
ax1.legend(['plasma','effect site','EC50'],fontsize=10)

ax2 = figure.add_subplot(2,1,2)
ax2.plot(times,effect, 'red',linewidth=2)
ax2.axis((0,480,0,1))
ax2.set_xlabel('time (min)', fontsize=12)
ax2.set_ylabel('effect (fraction of Emax)')
ax2.legend(['effect'],fontsize=10)

figure.savefig('standard_double_hyst_plot.png', dpi=300)  
figure.show-









    Out[46]:





<bound method Figure.show of <matplotlib.figure.Figure object at 0x0000000004040128>>

These plots use the following parameters:



In [50]:

    
t_half_effect = np.log(2)/ke0
(ke0, t_half_effect, ec50, nHill)









    Out[50]:





(0.00591003, 117.28319155062584, 336.28, 24.9464)

Units: 1/min, min, ng/ml, unitless



In [ ]: