ST 762, HOMEWORK 1, FALL 2009These problems are to be turned in on the due date.1. Many nonlinear (in parameters) functions used to describe biological and physical phenomenaarise as the solution to a system of ordinary differential equations (ODEs). This is the casein pharmacokinetics, as discussed in Chapter 1 of the notes, wh ere compartmental represen-tations of mechanisms taking place over time within a human or animal subject are routinelyused to characterize the pharmacokinetic processes of absorption, distribu tion, metabolism,and excretion. Similar representations are also used to describe processes underlying, forexample, the interplay between a virus such as human immunodeficiency virus (HIV) and theimmune system of an infected individual.On the bottom of page 22 of the class notes is a depiction of the one compartment modelwith first order absorption used to describe the pharmacokinetics of theophylline followingoral administration of dose D given at time 0 within a single individual. The compartmentalmodel may be translated into the following system of ODEs:dA(t)dt= kaAa(t) − keA(t)dAa(t)dt= −kaAa(t),(1)where A(t) is the amount of drug p resent in the main “blood” compartment at time t ≥ 0;and Aa(t) is the amount of drug present in a hypothetical “absorption depot” (e.g., the gut)at time t ≥ 0, from which it is absorbed into the main compartment at f ractional rate ofabsorption ka(units of 1/time). At time t = 0, it is assumed that the dose D instantaneouslyfills the “absorption depot,” so that we have the initial condition that Aa(t) = D. It is alsoassumed that there is no drug already pr esent in the system, so at time t = 0 there is nodrug in the main “blood” compartment; that is, we have th e initial condition A(0) = 0. Forsimplicity, we will take the b ioavailability F = 1.Because so many nonlinear mod els arise in this way, it is instructive to know something abouthow such s ystems of ODEs can be s olved (when they can be s olved in a closed form, thatis). Whether or n ot you have ever taken a cours e in differential equations, you will find bydoing a little research or by a little clever thinking that it is not too hard to solve this systemof equ ations under the given initial conditions to obtain the expr ession for the concentrationC(t) = A(t)/V present in the main compartment at time t given in (1.5) on page 23 of thenotes, where V is the hyp othetical volume of the main compartment. Your job in this problemis to give a full, step-by-step argument leading to an expression for the solution for A(t) usingany method you choose.Here are two possible ways to go about this (not the only ways):(a) Use the method of Laplace transforms, which you can research on line or in a standarddifferential equations text and apply directly to system (1).(b) Use “brute force” by carrying out the following steps: (i) solve the second equation in(1), which involves only Aa(t), for Aa(t) by integrating both sides of the equation andmaking use of the initial condition. (Hint: divide both sides of the equation by Aa(t)first.) (ii) Su bstitute the expression for Aa(t) into (i), rearrange the equation by placingall terms involving A(t) on the left hand side, multiply both sides by a su itable function1of t so that the new left hand side is the derivative of a product, integrate both sides,and make use of the initial condition.It is u p to you to choose a method of s olution. Present your argument leading to the solu-tion systematically and with suitable narrative so that a person unfamiliar with differentialequations could follow it.2. In your f avorite programming language, write two programs to implement the followin g meth-ods.(i) The 3-step GLS algorithm on pages 51 and 56 of the notes, where the “weight matrix”W is held fixed at step (iii), as d escribed in Section 3.2 of the notes. Your p rogramshould allow the user to choose the number of iterations C.(ii) IRWLS (C = ∞) as described in Section 3.4 of the notes.Do not mimic the programs in Section 3.7; I want you to write all parts of the algorithm (e.g.the Gauss-Newton scheme) yourself.Both programs should have the following features:• Allow the user to supply a starting value β(0)to get things going.• Allow the user to choose a maximum number of iterative updates in step (iii) of thefixed-C algorithm and similarly a maximum number of iterative updates of IRWLS; areasonable choice would be 500. If iteration continues up to this max number, eachprogram should stop and declare that no convergence was reached.• Use as the convergence criterion both in step (iii) of th e fixed-C algorithm and IRWLSthe following rule: If 2 successive iterates β(a+1)and β(a), s ay, have a relative differenceof less than some small constant tol, stop and declare β(a+1)to be the solution. That is,stop ifmaxℓ=1,...,p|βℓ,(a+1)− βℓ,(a)|/|βℓ,(a)| < tolas on page 59. For your programs, take tol = 10−8.• Compute an estimate of σ2based on the final estimated value for β. Use the versionof the estimator “adjusted for loss of degrees of freedom” (with the divisor (n − p)) onpage 64.To test your programs, consider the pharmacokinetic data fr om Subject 10 in the theophyllinestudy discussed in Example 1.8 in the class notes. Here, following oral dose of D = 5.50mg/kg, this s ubject had blood samples taken at times tj(hours), j = 1, . . . , 10, and on eachthe concentration of theophylline (mg/L), represented by the random variable Yjfor each j,was measured. Letting xj= (D, tj), the concentration-time relationship at this dose at timetjis thought to be well represented by the one compartment open model with fir st orderabsorption as in Problem 1, which we write in the formf(xj, β) =Deβ3eβ2(eβ3− eβ1/eβ2){exp(−eβ1tj/eβ2− exp(−eβ3tj)}, (2)Note that in (2), comparing to equation (1.5) of the class notes, th e model is parameter-ized in terms of β1= log Cl, β2= log V , and β3= log ka. The data are given in the file2theo10.dat, available on the class web page. The first column is time (hours) and the secondis concentration (mg/L).(a) Make a plot of the data using your favorite software. This will give you a sense of theshape of the concentration-time relationship for this subject.(b) Assume that the variance model isvar(Yj|xj) = σ2f2θ(xj, β), θ known.For this mean-variance model, try three different fits usin g each program:(i) θ = 0.0 (thus, assuming constant variance and fi