extantelementsextantelementsModeling Caffeine Metabolism - time offset<p>What’s the next level of complexity to add to our model? As noted previously, our model is not able to generate an arbitrary caffeine schedule. Let’s add functionality to enable it to do so!</p>
<p>What do we need in order to do this? What we would like to do is to be able to simulate a single day as a combination of caffeine inputs at arbitrary times during the day. We would then like to repeat the day periodically. To accomplish this, we need the ability to add a time offset to a dose delivered every 24 hours. This enables us to construct a day out of a combination of different doses with different time offsets, i.e. a 100 mg dose at 9 am and then a 50 mg dose at 4 pm.</p>
<p>Let’s first derive an expression for a periodic dose of <script type="math/tex">d</script> milligrams delivered every <script type="math/tex">a</script> hours with a time offset (delay) from <script type="math/tex">t=0</script> of <script type="math/tex">t_0</script>. As done previously, we can construct an expression for <script type="math/tex">f(t)</script>, the metabolically active amount at time <script type="math/tex">t</script>. We note that before <script type="math/tex">t_0</script>, the caffeine content will be zero as no dosing has started.</p>
<p>For <script type="math/tex">t \ge t_0</script>, this case is similar to the problem we have already solved but shifted in time by <script type="math/tex">t_0</script>. We can write out our sum of exponential terms (noting the time offset):</p>
<script type="math/tex; mode=display">f(t) = d\cdot e^{-(t-t_0)/\tau} + d\cdot e^{-(t-t_0 -a)/\tau} + d\cdot e^{-(t-t_0 - 2a)/\tau} + \ldots + d\cdot e^{-(t-t_0 - na)/\tau}</script>
<p>where <script type="math/tex">n</script> corresponds to the last consumption event before time <script type="math/tex">t</script>. We can figure out what <script type="math/tex">n</script> is by looking at the timing of events. There are no events before <script type="math/tex">t_0</script>. Afterwards, events occur every <script type="math/tex">a</script> hours. To account for the dead time caused by the offset, we simply subtract it from the total time and figure out how many events happened every <script type="math/tex">a</script> hours after that. Thus, we find that <script type="math/tex">n=\left\lfloor {\frac{t-t_0}{a}}\right\rfloor</script>.</p>
<p>We note that this expression is identical to our previous expression but with <script type="math/tex">t</script> replaced by <script type="math/tex">t-t_0</script>. Thus, the final sum will be:</p>
<script type="math/tex; mode=display">f(t\ge t_0) = d e^{-(t-t_0)/\tau} \left(\frac{1-e^{a(n+1)/\tau}}{1-e^{a/\tau}}\right)</script>
<p>Thus, we have our overall expression for <script type="math/tex">f(t)</script>:</p>
<script type="math/tex; mode=display">% <![CDATA[
f(t) = \left\{
\begin{array}{lr}
0 & t< t_0\\
d e^{-(t-t_0)/\tau} \left(\frac{1-e^{a(n+1)/\tau}}{1-e^{a/\tau}}\right) & t \geq t_0\\
\end{array}
\right. %]]></script>
<p>We have a simple expression for a repeated dose with a given offset. Now, we’d like the ability to build an arbitrary day out of such doses. So, let us consider doses <script type="math/tex">1, 2, \ldots, m</script> with dose size <script type="math/tex">d_1, d_2, \ldots, d_m</script>, offset <script type="math/tex">t_1, t_2, \ldots, t_m</script> and periodicity <script type="math/tex">a_1, a_2, \ldots, a_m</script>.</p>
<p>We can sum together these individual terms to construct the general expression. For <script type="math/tex">j=1,2,\ldots, m</script>:</p>
<script type="math/tex; mode=display">% <![CDATA[
f_j(t) = \left\{\begin{array}{lr}
0 & t<t_j\\
d_j e^{-(t-t_j)/\tau} \left(\frac{1-e^{a_j(n_j+1)/\tau}}{1-e^{a_j/\tau}}\right) & t \geq t_j\\
\end{array}\right. %]]></script>
<script type="math/tex; mode=display">f(t) = \sum_{j=1}^{m} f_j(t)</script>
<p>Note that the <script type="math/tex">n_j</script> values must be recalculated for every <script type="math/tex">j</script>.</p>
<p>To summarize so far, what we’ve done is constructed a sum of multiple different dosing schedules to be able to construct an arbitrary caffeine dosing schedule. We can implement these results in MATLAB as before. Note that we add a special case in the code for when <script type="math/tex">a=0</script>. This situation is nonphysical as it is equivalent to infinite consumption. We use <script type="math/tex">a=0</script> (durationBetween in the code) as a magic number to indicate a non-repeating dosage. We will return to this later. First, the code:</p>
<div class="highlight"><pre><code class="language-matlab" data-lang="matlab"><span class="k">function</span><span class="w"> </span>[ content ] <span class="p">=</span><span class="w"> </span><span class="nf">caffeineOffset</span><span class="p">(</span> time, caffeineContent, durationBetween, offset <span class="p">)</span><span class="w"></span>
<span class="w"> </span><span class="n">halfLife</span> <span class="p">=</span> <span class="mf">5.7</span><span class="p">;</span>
<span class="n">tau</span> <span class="p">=</span> <span class="n">halfLife</span><span class="o">/</span><span class="nb">log</span><span class="p">(</span><span class="mi">2</span><span class="p">);</span>
<span class="n">n</span> <span class="p">=</span> <span class="nb">floor</span><span class="p">((</span><span class="n">time</span><span class="o">-</span><span class="n">offset</span><span class="p">)</span><span class="o">/</span><span class="n">durationBetween</span><span class="p">);</span>
<span class="k">if</span> <span class="n">time</span> <span class="o"><</span> <span class="n">offset</span>
<span class="n">content</span> <span class="p">=</span> <span class="mi">0</span><span class="p">;</span>
<span class="k">elseif</span> <span class="n">durationBetween</span> <span class="o">==</span> <span class="mi">0</span>
<span class="n">content</span> <span class="p">=</span> <span class="n">caffeineContent</span> <span class="o">*</span> <span class="nb">exp</span><span class="p">(</span><span class="o">-</span><span class="p">(</span><span class="n">time</span><span class="o">-</span><span class="n">offset</span><span class="p">)</span><span class="o">/</span><span class="n">tau</span><span class="p">);</span>
<span class="k">else</span>
<span class="n">content</span> <span class="p">=</span> <span class="n">caffeineContent</span> <span class="o">*</span> <span class="nb">exp</span><span class="p">(</span><span class="o">-</span><span class="p">(</span><span class="n">time</span><span class="o">-</span><span class="n">offset</span><span class="p">)</span><span class="o">/</span><span class="n">tau</span><span class="p">)</span> <span class="o">*</span> <span class="p">((</span><span class="mi">1</span><span class="o">-</span><span class="nb">exp</span><span class="p">(</span><span class="n">durationBetween</span> <span class="o">*</span> <span class="p">(</span><span class="n">n</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span><span class="o">/</span><span class="n">tau</span><span class="p">))</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="nb">exp</span><span class="p">(</span><span class="n">durationBetween</span><span class="o">/</span><span class="n">tau</span><span class="p">)));</span>
<span class="k">end</span>
<span class="k">end</span></code></pre></div>
<div class="highlight"><pre><code class="language-matlab" data-lang="matlab"><span class="k">function</span><span class="w"> </span>[ content ] <span class="p">=</span><span class="w"> </span><span class="nf">caffeineOffsetList</span><span class="p">(</span> time, caffeineContent, durationBetween, offset <span class="p">)</span><span class="w"></span>
<span class="w"> </span><span class="c">%time is list of times</span>
<span class="c">%caffeineContent, durationBetween, and offset are lists of the same length </span>
<span class="n">content</span> <span class="p">=</span> <span class="nb">zeros</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span><span class="nb">length</span><span class="p">(</span><span class="n">time</span><span class="p">));</span>
<span class="k">for</span> <span class="nb">i</span> <span class="p">=</span> <span class="mi">1</span><span class="p">:</span><span class="nb">length</span><span class="p">(</span><span class="n">time</span><span class="p">)</span>
<span class="n">content</span><span class="p">(</span><span class="nb">i</span><span class="p">)</span> <span class="p">=</span> <span class="mi">0</span><span class="p">;</span>
<span class="k">for</span> <span class="nb">j</span> <span class="p">=</span> <span class="mi">1</span><span class="p">:</span><span class="nb">length</span><span class="p">(</span><span class="n">caffeineContent</span><span class="p">)</span>
<span class="n">content</span><span class="p">(</span><span class="nb">i</span><span class="p">)</span> <span class="p">=</span> <span class="n">content</span><span class="p">(</span><span class="nb">i</span><span class="p">)</span> <span class="o">+</span> <span class="n">caffeineOffset</span><span class="p">(</span><span class="n">time</span><span class="p">(</span><span class="nb">i</span><span class="p">),</span><span class="n">caffeineContent</span><span class="p">(</span><span class="nb">j</span><span class="p">),</span><span class="n">durationBetween</span><span class="p">(</span><span class="nb">j</span><span class="p">),</span><span class="n">offset</span><span class="p">(</span><span class="nb">j</span><span class="p">));</span>
<span class="k">end</span>
<span class="k">end</span>
<span class="k">end</span></code></pre></div>
<p>Let’s apply the code on some simple cases. The first thing we want to do is apply the model for a simple periodic dosing with some offset, just to confirm that the offset works correctly.</p>
<!--figure tags without plugin: http://stackoverflow.com/questions/19331362/using-an-image-caption-in-markdown-jekyll -->
<figure>
<a href="../images/Acaff_fig1.png">
<img src="../images/Acaff_fig1.png" alt="Figure 1: Response curve for a daily cup of coffee equivalent that is shifted 12 hours from t=0." />
</a>
<figcaption>
Figure 1: Response curve for a daily cup of coffee equivalent that is shifted 12 hours from t=0.
</figcaption>
</figure>
<p>It appears like the code works. Now we look at a slightly more interesting case, a sum of two different dosing schedule. This example has a 50 mg dose at midnight (<script type="math/tex">t=0,24,48,\ldots</script>) and a 100 mg dose at noon (<script type="math/tex">t=12,36,60,\ldots</script>), where we have set <script type="math/tex">t=0</script> as midnight. </p>
<!--figure tags without plugin: http://stackoverflow.com/questions/19331362/using-an-image-caption-in-markdown-jekyll -->
<figure>
<a href="../images/Acaff_fig2.png">
<img src="../images/Acaff_fig2.png" alt="Figure 2: Response curve for two doses a day at the given dosing pattern." />
</a>
<figcaption>
Figure 2: Response curve for two doses a day at the given dosing pattern.
</figcaption>
</figure>
<p>Now we can try something more complicated and crazy looking. The next example has a periodic 50 mg dose every 12 hours, with an initial offset of 18 hours, a daily 100 mg dose with an offset of 9 hours, and another daily 100 mg dose with an offset of 14 hours. This complicated dosing schedule is handled easily by the MATLAB code, resulting in Figure 3.</p>
<!--figure tags without plugin: http://stackoverflow.com/questions/19331362/using-an-image-caption-in-markdown-jekyll -->
<figure>
<a href="../images/Acaff_fig3.png">
<img src="../images/Acaff_fig3.png" alt="Figure 3: Response curve for the complex dosing schedule." />
</a>
<figcaption>
Figure 3: Response curve for the complex dosing schedule.
</figcaption>
</figure>
<p>The infinite time behavior is now more complex, as the individual periodic functions may interfere with each other due to varying <script type="math/tex">a_j</script> frequencies and <script type="math/tex">t_j</script> offsets.</p>
<p>The next thing we look at is the special feature mentioned just before the code–the ability to include non-repeating dosages. This allows the program to act as a logbook of sorts to see what one’s caffeine levels might be. In this case, a list of inputs and offsets are provided where the offsets span the total duration of the time period of interest. In this case, we are no longer summing over many terms, so we can remove the summed geometric series from the expression and just replace it with the simpler exponential term.</p>
<p>As an example, the doses <script type="math/tex">[100,100,50,50,100,50]</script> are associated with the offsets <script type="math/tex">[9,15,32,45,65,80]</script>, respectively. The resulting curve is plotted in Figure 4.</p>
<!--figure tags without plugin: http://stackoverflow.com/questions/19331362/using-an-image-caption-in-markdown-jekyll -->
<figure>
<a href="../images/Acaff_fig4.png">
<img src="../images/Acaff_fig4.png" alt="Figure 4: Logbook like behavior." />
</a>
<figcaption>
Figure 4: Logbook like behavior.
</figcaption>
</figure>
<p>To summarize the results of this post: we expanded our method to handle arbitrary offsets and allowed for superposition of multiple periodic dosing regimens (and single doses) in order to be able to construct any possible dosing scheme. We looked at several different examples of periodic dosing schemes as well as completely manual logbook-style calculations.</p>
<p><strong>Again, it must be stressed again that this is a toy model. Real systems don’t work this way, due to many physical effects and feedback systems within a person. It should also be noted that caffeine is used throughout as an example of interest, but it could be any other pharmaceutical drug that follows first-order exponential decay kinetics.</strong></p>
Fri, 09 Jan 2015 00:00:00 +0000
/Modeling-Caffeine-Metabolism---Time-Offset/
/Modeling-Caffeine-Metabolism---Time-Offset/Modeling Caffeine Metabolism - mathematical note<p>I wanted to make a quick mathematical note on the long term behavior of our caffeine metabolism function. In the previous post, we came up with the expressions </p>
<p><script type="math/tex">f(t)=d\cdot e^{-t/\tau}</script> as the exponential decay for a single dose of <script type="math/tex">d~</script> milligrams at time <script type="math/tex">t</script>.</p>
<p><script type="math/tex">\displaystyle f(t) = d e^{-t/\tau} \left(\frac{1-e^{a(n+1)/\tau}}{1-e^{a/\tau}}\right)</script> for the behavior at time <script type="math/tex">t~</script> for a dose of <script type="math/tex">d~</script> mg delivered every <script type="math/tex">a~</script> hours, where <script type="math/tex">n=\lfloor t/a \rfloor</script>.</p>
<p>In this post, we’ll look at the long term behavior using two similar approaches. The quantity we are most interested in is the amount of active caffeine that is always present regardless of where the person is along the dosing cycle, i.e. the lowest point along the exponential curve (just before the spike up).</p>
<p><strong>Approach 1 - Infinite Series</strong></p>
<p>We know that just before <script type="math/tex">t=a</script> (just before the spike), we have that the amount is given by <script type="math/tex">d e^{-a/\tau}</script>. Similarly, just before <script type="math/tex">t=2a</script>, we have that the amount is given by <script type="math/tex">de^{-a/\tau} + de^{-2a/\tau}</script>. Generalizing to <script type="math/tex">t=\infty a</script>, we have the amount given by the infinite series</p>
<script type="math/tex; mode=display">\displaystyle T = de^{-a/\tau} + de^{-2a/\tau} + de^{-3a/\tau} + \ldots = \frac{de^{-a/\tau}}{1-e^{-a/\tau}}</script>
<p>where we have summed the infinite series. This expression gives us the desired quantity.</p>
<p><strong>Approach 2 - Extending Finite Series</strong></p>
<p>We are interested in the behavior at times just before a new dose. We can represent these times as <script type="math/tex">t=ma-\epsilon</script> where <script type="math/tex">m</script> is a positive integer and <script type="math/tex">\epsilon \ll a</script>. We use <script type="math/tex">\epsilon</script> so that we are working just before the new dose, occurring at <script type="math/tex">t=ma</script>.</p>
<p>We can also calculate <script type="math/tex">n=\left\lfloor \frac{ma-\epsilon}{a}\right\rfloor = \left\lfloor m - \frac{\epsilon}{a}\right\rfloor = m-1</script>.</p>
<p>We can substitute into our formula for <script type="math/tex">f(t)</script> .</p>
<script type="math/tex; mode=display">f(ma-\epsilon) = d e^{-(ma-\epsilon)/\tau} \left(\frac{1-e^{am/\tau}}{1-e^{a/\tau}}\right)</script>
<script type="math/tex; mode=display">f(ma-\epsilon) = d \left(\frac{e^{-ma/\tau} - e^{-am/\tau}e^{am/\tau}}{1-e^{a/\tau}} \right) e^{\epsilon/\tau}</script>
<script type="math/tex; mode=display">f(ma-\epsilon) = d\left( \frac{e^{-ma/\tau} - 1}{1-e^{a/\tau}} \right) e^{\epsilon/\tau}</script>
<script type="math/tex; mode=display">f(ma-\epsilon) = d\left(\frac{e^{-ma/\tau} - 1}{e^{a/\tau}(e^{-a/\tau}-1)} \right) e^{\epsilon/\tau}</script>
<script type="math/tex; mode=display">f(ma-\epsilon) = de^{-a/\tau}\left(\frac{1-e^{-ma/\tau}}{1-e^{-a/\tau}} \right) e^{\epsilon/\tau}</script>
<p>Now we consider the limit we are working in–unbounded <script type="math/tex">m</script> and small <script type="math/tex">\epsilon</script>. As <script type="math/tex">m\rightarrow \infty</script>, <script type="math/tex">e^{-ma/\tau} \rightarrow 0</script> and as <script type="math/tex">\epsilon\rightarrow 0</script>, <script type="math/tex">e^{\epsilon/\tau} \rightarrow 1</script>. Accordingly, we have this expression approaching <script type="math/tex">\displaystyle\frac{de^{-a/\tau}}{1-e^{-a/\tau}}</script> as found previously.</p>
<p>Through two slightly different approaches, we were able to determine the long term behavior of the function. For very large <script type="math/tex">m</script>, we have that the function spikes up to</p>
<p><script type="math/tex">\displaystyle \frac{d}{1-e^{-a/\tau}} </script> at <script type="math/tex">t=ma</script> </p>
<p>and decays down to</p>
<p><script type="math/tex">\displaystyle \frac{de^{-a/\tau}}{1-e^{-a/\tau}}</script> just before <script type="math/tex">t=(m+1)a</script>.</p>
Sat, 03 Jan 2015 00:00:00 +0000
/Modeling-Caffeine-Metabolism--mathematical-note/
/Modeling-Caffeine-Metabolism--mathematical-note/Modeling Caffeine Metabolism<p>Caffeine is a stimulant drug that is widely consumed in tea, coffee, soda, and chocolate, among other beverages and foods. In this post, we construct a toy model for caffeine metabolism.</p>
<p>The mechanism of action of caffeine is complex, as is the process of developing tolerance. These are important factors to consider when thinking about the ultimate effect of the drug on alertness and we will discuss them briefly at the end.</p>
<p>So–let’s build our model. We need to consider caffeine input (consumption) and caffeine output (metabolism). </p>
<p>For our caffeine input, let’s consider a semi-realistic model of consuming <script type="math/tex">d</script> milligrams of caffeine as a single dose every <script type="math/tex">a</script> hours. The simplifying assumptions associated with this model include the non-physical instant consumption of an entire source of caffeine under a rigid timeframe as well as neglecting the time needed for absorption within the body.</p>
<p>The instant consumption is not a terrible approximation since most sources are often consumed within 15 minutes or so. The rigid timefame is necessary to make the mathematics simple (as will be shown later), and is a more significant restriction. This model could not be used to consider someone who consumes a cup of coffee only at 9:00 am and 3:00pm every day, due to the periodic constraint on the caffeine input. This restriction imposes a constraint on the problems we can look at.</p>
<p>The assumption on the absorption of caffeine is a significant approximation. Typically, caffeine is absorbed through the intestine within 45 minutes of consumption with the peak blood caffeine content occuring between 1-2 hours after consumption. We have ignored these effects governing the bioavailability of caffeine. To justify this assumption, we can say that what our model is really looking at is the blood caffeine content with an instant input, i.e. instanteous absorption from the intestine into bioavailable forms that is time-shifted by some consistent amount after actual consumption. However, in a real system, there will be non-trivial dispersion of this input.</p>
<p>For our caffeine output (metabolism), we model the remaining active amount of the drug as an exponential decay with half-life 5.7 hours, based on the literature cited below. These assumptions are reasonable, given experimental data on caffeine content. Ignoring any prefactors, the significant tunable parameter for an exponential decay curve is the half-life (or equivalently, the time-constant). The half-life of caffeine in humans has been demonstrated to vary wildly depending on a variety of states of health. 5.7 hours is a reasonable mean value. The biochemistry underlying this decay is the degradation of caffeine into a number of other biomolecules, and follows first order kinetics.</p>
<p>We can now begin to put together our mathematical model. </p>
<p>First, the half-life of caffeine can be related to the time constant (for convenience) by</p>
<script type="math/tex; mode=display">\lambda = \tau \cdot \ln{2}</script>
<script type="math/tex; mode=display">\tau=\frac{\lambda}{\ln{2}} = \frac{5.7}{\ln{2}}\textrm{hrs.} = 8.223 \textrm{ hrs.}</script>
<p>This gives our simple exponential decay for a single dose of <script type="math/tex">d</script> mg at time <script type="math/tex">t=0</script> as <script type="math/tex">f(t)=d\cdot e^{-t/\tau}</script> where <script type="math/tex">f(t)</script> gives the amount of metabolically active caffeine.</p>
<p>We can now evaluate a sum of these terms that are each time-shifted by the time of intake of caffeine. As a simple example, for caffeine consumption at <script type="math/tex">t=0</script> and <script type="math/tex">t=a</script>, we have</p>
<script type="math/tex; mode=display">f(t) = d\cdot e^{-t/\tau} + d\cdot e^{-(t-a)/\tau}</script>
<p>for all time <script type="math/tex">t\geq a</script> since the second dose must be consumed for the equation to be valid.</p>
<p>For a more general sum at time <script type="math/tex">t=t_0</script>, we need to sum the terms for all consumption events that have happened before <script type="math/tex">t_0</script>. These events occur at time <script type="math/tex">t=0,a,2a,\ldots, n\cdot a</script> where <script type="math/tex">na</script> is the largest multiple of <script type="math/tex">a</script> that is less than or equal to <script type="math/tex">t_0</script>. We find that <script type="math/tex">n=\left\lfloor {\frac{t_0}{a}}\right\rfloor</script>. Thus, we can write our sum as a sum of these decaying exponentials, with the number of exponential terms equal to the number of consumption events that have occured:</p>
<script type="math/tex; mode=display">f(t) = \sum_{i=0}^{n}{d\cdot e^{-(t-i\cdot a)/\tau}}</script>
<script type="math/tex; mode=display">f(t) = d \sum_{i=0}^{n}{e^{-t/\tau} e^{i\cdot a/\tau}}</script>
<script type="math/tex; mode=display">f(t) = d e^{-t/\tau} \sum_{i=0}^{n}{e^{i\cdot a/\tau}}</script>
<p>We recognize the term within the summation to be a geometric series and replace the summation with the closed form sum.</p>
<script type="math/tex; mode=display">f(t) = d e^{-t/\tau} \left(\frac{1-e^{a(n+1)/\tau}}{1-e^{a/\tau}}\right)</script>
<p>We now have our closed form expression for the active caffeine content. A direct observation that we can make is that the consumption dose acts as a scalar multiplier over the entire function.</p>
<p>What can we do with this expression? Well, we can use it to model several profiles of caffeine consumption. I implemented the model in MATLAB. The <strong>caffeine</strong> function takes in a single time point, the caffeine content of the source <script type="math/tex">d</script>, and the duration between consumptions events <script type="math/tex">a</script> and returns the active caffeine amount at that time. </p>
<div class="highlight"><pre><code class="language-matlab" data-lang="matlab"><span class="k">function</span><span class="w"> </span>[ content ] <span class="p">=</span><span class="w"> </span><span class="nf">caffeine</span><span class="p">(</span> time, caffeineContent, durationBetween <span class="p">)</span><span class="w"></span>
<span class="w"> </span><span class="n">halfLife</span> <span class="p">=</span> <span class="mf">5.7</span><span class="p">;</span>
<span class="n">tau</span> <span class="p">=</span> <span class="n">halfLife</span><span class="o">/</span><span class="nb">log</span><span class="p">(</span><span class="mi">2</span><span class="p">);</span>
<span class="n">n</span> <span class="p">=</span> <span class="nb">floor</span><span class="p">(</span><span class="n">time</span><span class="o">/</span><span class="n">durationBetween</span><span class="p">);</span>
<span class="n">content</span> <span class="p">=</span> <span class="n">caffeineContent</span> <span class="o">*</span> <span class="nb">exp</span><span class="p">(</span><span class="o">-</span><span class="n">time</span><span class="o">/</span><span class="n">tau</span><span class="p">)</span> <span class="o">*</span> <span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="nb">exp</span><span class="p">(</span><span class="n">durationBetween</span> <span class="o">*</span> <span class="p">(</span><span class="n">n</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span><span class="o">/</span><span class="n">tau</span><span class="p">))</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="nb">exp</span><span class="p">(</span><span class="n">durationBetween</span><span class="o">/</span><span class="n">tau</span><span class="p">));</span>
<span class="k">end</span></code></pre></div>
<p>The <strong>caffeineList</strong> function takes in a list of times, <script type="math/tex">d</script>, and <script type="math/tex">a</script>, and returns a list containing the active caffeine amounts at the times indicated by the input.</p>
<div class="highlight"><pre><code class="language-matlab" data-lang="matlab"><span class="k">function</span><span class="w"> </span>[ content ] <span class="p">=</span><span class="w"> </span><span class="nf">caffeineList</span><span class="p">(</span> time, caffeineContent, durationBetween <span class="p">)</span><span class="w"></span>
<span class="w"> </span><span class="n">content</span> <span class="p">=</span> <span class="nb">zeros</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span><span class="nb">length</span><span class="p">(</span><span class="n">time</span><span class="p">));</span>
<span class="k">for</span> <span class="nb">i</span> <span class="p">=</span> <span class="mi">1</span><span class="p">:</span><span class="nb">length</span><span class="p">(</span><span class="n">time</span><span class="p">)</span>
<span class="n">content</span><span class="p">(</span><span class="nb">i</span><span class="p">)</span> <span class="p">=</span> <span class="n">caffeine</span><span class="p">(</span><span class="n">time</span><span class="p">(</span><span class="nb">i</span><span class="p">),</span><span class="n">caffeineContent</span><span class="p">,</span><span class="n">durationBetween</span><span class="p">);</span>
<span class="k">end</span>
<span class="k">end</span></code></pre></div>
<p>The first thing we should do with this model is the trivial case of simulating a single dose. This acts as a sanity check to make sure that things look reasonable. This yields:</p>
<!--figure tags without plugin: http://stackoverflow.com/questions/19331362/using-an-image-caption-in-markdown-jekyll -->
<figure>
<a href="../images/caff_fig1.png">
<img src="../images/caff_fig1.png" alt="Figure 1: Response curve for a single cup of coffee equivalent at time 0." />
</a>
<figcaption>
Figure 1: Response curve for a single cup of coffee equivalent at time 0.
</figcaption>
</figure>
<p>Cool! It looks like the model works on this trivial case. Now we can look at more interesting cases. The next example is a coffee-cup-equivalent every 24 hours.</p>
<!--figure tags without plugin: http://stackoverflow.com/questions/19331362/using-an-image-caption-in-markdown-jekyll -->
<figure>
<a href="../images/caff_fig2.png">
<img src="../images/caff_fig2.png" alt="Figure 2: Response curve for a daily cup of coffee equivalent." />
</a>
<figcaption>
Figure 2: Response curve for a daily cup of coffee equivalent.
</figcaption>
</figure>
<p>We can now look at someone who drinks coffee twice a day. Unfortunately, we’re forced by our formalism to stick to doses that are 12 hours apart, which is perhaps not the most realistic model for how people consume coffee, but still perhaps worthwhile to look at.</p>
<!--figure tags without plugin: http://stackoverflow.com/questions/19331362/using-an-image-caption-in-markdown-jekyll -->
<figure>
<a href="../images/caff_fig3.png">
<img src="../images/caff_fig3.png" alt="Figure 3: Response curve for two cups of coffee equivalent per day." />
</a>
<figcaption>
Figure 3: Response curve for two cups of coffee equivalent per day.
</figcaption>
</figure>
<p>Now, we look at frequent consumption of a smaller dose of caffeine. Our model here is four servings of 54 mg of caffeine, simulating the extreme consumption of four cans of soda per day.</p>
<!--figure tags without plugin: http://stackoverflow.com/questions/19331362/using-an-image-caption-in-markdown-jekyll -->
<figure>
<a href="../images/caff_fig4.png">
<img src="../images/caff_fig4.png" alt="Figure 4: Reponse curve for four cans of soda equivalent per day." />
</a>
<figcaption>
Figure 4: Reponse curve for four cans of soda equivalent per day.
</figcaption>
</figure>
<p>Now we can think about the results a bit. First, the shape of the curves is entirely unsurprising–this is what a train of exponentials looks like. The result may be familiar from circuit theory. </p>
<p>We have two time-scales to look at: the steady-state and the transient. </p>
<p>The steady-state is perhaps more interesting, simulating the behavior of a habitutated drinker. We can see that in the case of 1 cup of coffee per day (figure 2), the caffeine content actually dips quite low before the next dose. However, when the drinker goes up to two cups per day (figure 3), they have a significant background level of active caffeine that is around 60 mg. This means that that much caffeine is constantly active in their bodies. </p>
<p>In the case of the transients, we can think of this as someone coming back to caffeine consumption after a long period without any consumption. The transition from no consumption to steady-state is perhaps clearest in figure 4. Here, we see the rise in background up to about 50 mg of always active caffeine (nearly a can of soda equivalent).</p>
<p>Remember, this is a toy model and shouldn’t be used for quantitative predictions. Qualitatively, we can clearly see that frequent consumption leads to a significant background level.</p>
<p>We can briefly consider the biological side of things. The mechanism of action of caffeine is to act as a competitive inhibitor on certain receptors within certain neurons. An increased background level of cafffeine means that the caffeine is constantly acting on these receptors as a competitive inhibitor. In a realistic biological system, this activity may change the profile of gene expression to increase receptor production, requiring more caffeine to result in the same effect, i.e. tolerance to a certain dosage of caffeine.</p>
<p>It is important to note that the relationship between our computed curves and complex physiological effects such as increased alertness is very complex. From our model, we really cannot say much with any certainty about this relationship, due to the complexity of the relationships between amounts of the drug, its biochemical activity, and its effect on behavior. There may be complex feedback regulation at multiple levels that makes this relationship complex and beyond the scope of this simple model.</p>
<p>To summarize: we have built a toy model for caffeine metabolism and used it to evaluate several different profiles of caffeine consumption, with a perspective towards the underlying biochemistry.</p>
<p><strong>Disclaimer: this is not medical advice. The information in this post is not intended to diagnose, treat, cure, or prevent any disease.</strong></p>
<p><strong>References</strong></p>
<ol>
<li>“Serum caffeine half-lives. Healthy subjects vs. patients having alcoholic hepatic disease.” Am J Clin Pathol. 1980 Mar;73(3):390-3. <a href="http://www.ncbi.nlm.nih.gov/pubmed/7361718">PubMed</a></li>
<li>“Plasma and salivary pharmacokinetics of caffeine in man.” Eur J Clin Pharmacol. 1981;21(1):45-52. <a href="http://www.ncbi.nlm.nih.gov/pubmed/7333346">PubMed</a></li>
<li>Wikipedia page on <a href="https://en.wikipedia.org/wiki/Caffeine">Caffeine</a>.</li>
</ol>
Tue, 30 Dec 2014 00:00:00 +0000
/Modeling-Caffeine-Metabolism/
/Modeling-Caffeine-Metabolism/this is a test of math and code<script type="math/tex; mode=display">i\hbar \frac{\partial}{\partial t}\psi=\hat{H}\psi</script>
<div class="highlight"><pre><code class="language-c" data-lang="c"><span class="kt">int</span> <span class="nf">main</span><span class="p">(</span> <span class="p">)</span> <span class="p">{</span>
<span class="n">printf</span><span class="p">(</span><span class="s">"hello, world"</span><span class="p">);</span>
<span class="k">return</span> <span class="mi">0</span><span class="p">;</span>
<span class="p">}</span></code></pre></div>
Tue, 30 Dec 2014 00:00:00 +0000
/Math-and-Code-Test/
/Math-and-Code-Test/hello, world<p>hello, world</p>
Tue, 30 Dec 2014 00:00:00 +0000
/Hello-World/
/Hello-World/