An adaptive drug delivery design using neural networks for effective treatment of infectious diseases: A simulation study

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 4 ( 2 0 0 9 ) 207–222

journal homepage: www.intl.elsevierhealth.com/journals/cmpb

An adaptive drug delivery design using neural networks for effective treatment of infectious diseases: A simulation study Radhakant Padhi a,∗ , Jayender R. Bhardhwaj b a b

Department of Aerospace Engineering, Indian Institute of Science-Bangalore, Karnataka 560012, India Department of Aerospace Engineering, Indian Institute of Technology-Madras, 600036, India

a r t i c l e

i n f o

a b s t r a c t

Article history:

An adaptive drug delivery design is presented in this paper using neural networks for effec-

Received 15 October 2007

tive treatment of infectious diseases. The generic mathematical model used describes the

Received in revised form

coupled evolution of concentration of pathogens, plasma cells, antibodies and a numeri-

9 December 2008

cal value that indicates the relative characteristic of a damaged organ due to the disease

Accepted 9 December 2008

under the influence of external drugs. From a system theoretic point of view, the external drugs can be interpreted as control inputs, which can be designed based on control theo-

Keywords:

retic concepts. In this study, assuming a set of nominal parameters in the mathematical

Infectious disease

model, first a nonlinear controller (drug administration) is designed based on the principle

Immunology dynamics

of dynamic inversion. This nominal drug administration plan was found to be effective in

Drug delivery

curing “nominal model patients” (patients whose immunological dynamics conform to the

Dynamic inversion

mathematical model used for the control design exactly. However, it was found to be inef-

Neuro-adaptive design

fective in curing “realistic model patients” (patients whose immunological dynamics may

Biomedical control

have off-nominal parameter values and possibly unwanted inputs) in general. Hence, to make the drug delivery dosage design more effective for realistic model patients, a modelfollowing adaptive control design is carried out next by taking the help of neural networks, that are trained online. Simulation studies indicate that the adaptive controller proposed in this paper holds promise in killing the invading pathogens and healing the damaged organ even in the presence of parameter uncertainties and continued pathogen attack. Note that the computational requirements for computing the control are very minimal and all associated computations (including the training of neural networks) can be carried out online. However it assumes that the required diagnosis process can be carried out at a sufficient faster rate so that all the states are available for control computation. © 2009 Elsevier Ireland Ltd. All rights reserved.

1.

Introduction

The immune system of living organisms exists to defend it from agents with properties of genetically alien information (such as bacteria, viruses, proteins, tissue and transformed

∗

own cells such as tumor cells). For this purpose, there exist three distinct levels of defense against the invading microbes. The first level is the outer perimeter of defense, which is the surface epithelial layers of the body. This includes the epidermal cells of the skin and the mucosal cells that line the

Corresponding author. Tel.: +91 80 2293 2756; fax: +91 80 2360 0134. E-mail address: [email protected] (R. Padhi). 0169-2607/$ – see front matter © 2009 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.cmpb.2008.12.010

208

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 4 ( 2 0 0 9 ) 207–222

respiratory, gastrointestinal and genitourinary tracts. Once this type of defence is surpassed, the innate immune system provides a tactical response, signaling the presence of non-self organisms and activating B-cells to produce antibodies that bind to the intruder’s antigens, reducing them to non-functioning units. They also stimulate the production of molecules that either damage the intruder’s plasma membrane directly or help trigger the second phase of immune response. The innate immune system protects against many extra cellular bacteria or free viruses that are found in blood plasma, lymph, tissue fluid, or interstitial space between cells. The adaptive immune system produces protective cells that remember specific antigens and produce antibodies customized to the pathogens. The generation of such defending cells and molecules by the immune system is termed as the immune response. An interested reader can refer to [2,5,12,21,24,25,28,29] for better understanding of the associated complex mechanism. Whenever there is an invasion of pathogens in an organ, the immune system responds through the production of plasma cells. The plasma cells produce antibodies that fight the pathogen. If the concentration of the invading pathogen is high, there is considerable damage to the organ before the body can launch a strong immune response. As a result, the immune system becomes weak and it leads to a lowered production of plasma cells and antibodies. When this happens, the affected organ fails to recover naturally and there is a need for external medication. The options available for treatment of an infectious disease include killing the invading microbes directly, enhancing the efficacy of the immune system (i.e. production of plasma cells and antibodies) or providing healing care to the affected organ directly [34–36]. Such drugs generally augment the natural ability of the body to fight the infection and eventually cure the damaged organ. However, in this work we have assumed only the availability of drugs that kill the invading microbes and heal the affected organ. We have not considered the drugs that enhance the efficacy of the immune system with the the following motivations: (i) the drugs that enhance the production of plasma cells and antibodies are not readily available (it is still a topic of pharmaceutical research), (ii) the drugs we have considered are effective enough to achieve our objectives, (iii) elimination of extra drugs also leads to the elimination of sideeffects and (iv) the treatment plan possibly becomes more cost effective. Unfortunately, all drugs come with the price of unwanted side-effects. Moreover, the response of the body to any external drug is not instantaneous. Hence a carefully-monitored continuous low dosage is preferable as compared to high dosage at discrete points in time (impulse control), which is currently in practice. However, since it is impossible for a doctor to monitor a patient and decide an appropriate quantity of the drug in a continuous manner, the idea can be realized in practice only if the process is made automatic (monitored and controlled by a computer). This is where control systems theory becomes a valuable tool. Using the various advanced control design methods, an appropriate control (drug) requirement history can be found out, which is both effective in treating the disease and distributed over time, leading to minimal side effects.

Dynamic models of the immune system have been described as early as 1978 by Perelson [5], where a set of differential equations have been developed to model the dynamics of antibodies in mammals. The concept of binary strings was used to describe and measure the strength of reactions between antibodies and antigens. In their approach, the number of state variables is also a function of time, taking into account the death of antibodies due to self-destruction and formation of new types of antibodies. In 1994, Asachenkov et al. [1] described the human immunological dynamics using a generic nonlinear mathematical model. Based on the idea presented in [1], a model including the effect of various drug (control) inputs is presented in [34–36], which is then used to synthesize a meaningful controller based on well-established classical optimal control theory. In [34], a nonlinear optimal controller is obtained by solving the associated two-point boundary value problem, using the “gradient method” [16]. This work is followed up in [35], where the authors have augmented this controller with a linearized “neighboring optimal controller”, which is based on standard linear quadratic regulator (LQR) theory and come up with a state feedback design for the linearized deviation dynamics. This work is further followed up in [36], where the standard linear quadratic gaussian (LQG) theory is used to address the issue of process and measurement noise. LQG design essentially combines the LQR theory and Kalman filter design to estimate the unmeasured and/or corrupted states. Compared to optimal control design, a relatively simpler and popular method of nonlinear control design is the technique of “dynamic inversion”, which is essentially based on the philosophy of feedback linearization [14,32]. In this approach, an appropriate co-ordinate transformation is carried out to make the system dynamics take a linear form in the transformed coordinates. Then linear control design tools are used to synthesize the controller. This technique guarantees global asymptotic (rather exponential) stability of the output error dynamics. Based on this dynamic inversion theory, in this study we propose an alternate drug delivery design approach, which has the following characteristics: (i) it comes up with a state feedback closed form solution. This makes the approach computationally non-intensive, and hence it can be implemented in real time easily; (ii) the control strategy is made adaptive to the actual parameter values of individual cases, thereby tailoring the drug delivery and making it patient-specific and (iii) there is no need to predict the treatment duration apriori (i.e. the medication can be continued for arbitrarily long period of time and can be stopped whenever the condition of the patient is found to improve. A major drawback of the dynamic inversion approach, however, is its sensitivity to modeling errors and parameter inaccuracies [4]. To address this issue, Calise and his co-workers have proposed to augment the dynamic inversion technique with neural networks so that the inversion error is canceled out (see [3,9,15–18] for example). This augmented technique makes the dynamic inversion technique more robust to modeling and parameter inaccuracies and hence makes it an extremely useful technique. Despite its significant contribution and popularity, the basic approach followed by Calise and his co-workers mainly serves as an augmenting tool for the dynamic inversion design, and hence it

209

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 4 ( 2 0 0 9 ) 207–222

cannot be used in conjunction with other promising control design approaches (optimal control design, for example). In contrast, in this study we have followed a model-following neuro-adaptive approach which can be used in conjuction with “any” established control design technique (e.g. dynamic inversion, optimal control, back-stepping, etc.). Note that the neuro-adaptive approach followed here can also be thought of as an augmentation tool to the approaches followed in [34–36], to make the control design scheme proposed there robust to the issues of parametric inaccuracy and modeling errors. There has been a phenomenal growth in the field of artificial neural networks and its applications to control systems in the last three decades. A compiled book is available detailing various applications of artificial neural networks with respect to identification and control of dynamic systems [23]. A good survey paper [11] outlines various applications of neural networks to control system design. An interested reader can refer to these references and references therein for a comprehensive knowledge about some of the early developments in the field. One potential application for using neural networks in control design is to identify and compensate for plant nonlinearities. In [19] such an approach is presented taking the help of multi-layer feed-forward neural networks that leads to a better tracking performance. A direct adaptive tracking control architecture with Gaussian radial basis function networks to compensate for plant nonlinearities is presented in [30], where the authors showed that neural networks can infact be trained online, and hence, require no prior training. The update process proposed in [30] also keeps the weights of the neural networks bounded. In a recent paper [31], an adaptive output feedback control scheme for the output tracking of a class of nonlinear systems using radial basis function (RBF) neural networks is presented, where the weights are adapted using a Lyapunov-based design as well. A systematic model-following neuro-adaptive approach appeared in [10] for a class of nonlinear systems, which comes up with an “extra control”, which when added to any nominal controller makes it robust to parameter uncertainties and a class of modeling inaccuracies. In a very recent paper [26], an indirect model-following neuro-adaptive control design approach is presented for addressing the same issues for non-square, non-affine systems. The method used here relies on this idea. In this study, first we have used the dynamic inversion technique to design a nominal controller (i.e. the control which is designed assuming that the model has nominal parameter values). We have also proved that the associated internal dynamics [14,32] remains stable with the application of this control. Next, the control strategy is made adaptive to the actual parameter values of individual cases with the help of neural networks trained online, thereby tailoring the drug delivery and making it patient-specific. Numerical results clearly demonstrate that the resulting drug delivery scheme (controller) is adaptive both with respect to the parameter uncertainties and also with respect to an unknown rate of continuing introduction of pathogens. Rest of the paper is organized as follows. In Section 2, the mathematical model for the immunological dynamics used in this work is described, which is followed the goal for the proposed drug delivery design. The system theory based mathematical details are presented in Section 3. Extensive

simulation results are contained in Section 4, which show that the proposed adaptive automatic drug delivery strategy can leads to complete recovery of a realistic model patient, even with large parameter inaccuracies and/or continuing introduction of pathogens. Appropriate conclusions are drawn in Section 5.

2. Immunology dynamics and objective of the control design 2.1.

Immunology dynamics

A generic nonlinear mathematical model describing the human immunology dynamics is presented in [1]. The model describes the coupled evolution of concentration of pathogens, plasma cells and antibodies and a numerical value that indicates the relative characteristic of a damaged organ. Stengel et al. [34–36] have augmented this model to incorporate the effect of four generic drugs (control inputs); namely pathogen killer, plasma cell enhancer, antibody enhancer and organ healing factor. An interested reader can see these references for more details. However, in this study we have assumed only the availability of drugs that kill the invading microbes and heal the affected organ. We have not considered drugs that enhance the efficacy of the immune system for the reasons pointed out in Section 1. With the effects of these drugs, the nonlinear mathematical model for the immunology dynamics can be written as ˙ d = fd (Xd ) + B Ud X

(1)

where Xd = [x1d x2d x3d x4d ]T , Ud = [u1d u2d ]T are the state and control vectors respectively. Function fd (Xd ) and matrix B are defined as

⎡ fd (Xd ) 

⎤

(a11d − a12d x1d )x3d ⎢ a21d (x4d )a22d x1d x3d − a23d (x2d − x2∗ ) ⎥

⎢ ⎣

⎡

b1 ⎢0 B ⎣ 0 0

a31d x2d − (a32d + a33d x1d )x3d a41d x1d − a42d x4d

d

⎥, ⎦

⎤

0 0 ⎥ ⎦ 0 b4

(2)

Here x1d represents concentration of the pathogen, x2d represents concentration of plasma cells (which are carriers and producers of antibodies), x3d represents concentration of antibodies (which kill the pathogen) and x4d represents relative characteristic of a damaged organ (or organ health factor). Similarly u1d represents the drug which is a pathogen killer and u2d represents the drug which heals the affected organ. x2∗ represents the steady state concentration of plasma cells. d The parameter a21d is a nonlinear function of x4d that describes the immune deficiency caused by damage to the organ and is given by

 a21d (x4d ) =

cos(x4d ), 0 ≤ x4d ≤ 1/2 0, x4d > 1/2

(3)

210

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 4 ( 2 0 0 9 ) 207–222

Table 1 – Nominal parameter values. Parameter

Value

Parameter

Value

a11d a12d a21d a22d a23d a31d

1 1 See (3) 3 1 1

a32d a33d a41d a42d b1 b4

1.5 0.5 0.5 1 −1 −1

The reason for introducing the subscript ‘d’ arises from the fact that it represents the “desired dynamics” in our formulation. The necessity of this will be clear later in this section. The numerical values of other parameters are as given in Table 1[34–36]. In the model, the parameter a11d indicates exponential growth of the pathogens. The term a12d x3d x1d designates the number of antigens neutralized by the antibodies x3d , where a12d is the coefficient related to the probability of neutralization of the germs by the antibodies upon an encounter. The term a41d x1d represents the degree of damage to the organ by the pathogens while the term a42d x4d signifies the recuperative capacity of the organ. For more information about the physical meanings of the various terms of the model, the reader is referred to [1,22,34–36]. Note that the model in (1)–(3)(with the associated numerical values of the parameters in Table 1) is a generic one and the numerical values in the model ‘qualitatively’ describe the various effects observed in practice in general. The main motivation of this study, however, is to present a general nonlinear adaptive control design technique for automatic treatment of infectious diseases. The model and parameter values considered are sufficient to demonstrate this. For specific diseases the model will probably remain the same, but it will have different parameter values. However, it is important to note that as and when such a model is available (to the best of the knowledge of authors such a model is not available yet), the nonlinear control design technique presented in this paper can directly be applied on it. To demonstrate the effectiveness of the adaptive control design strategy presented in this paper, however, we have assumed that the generic model in (1)–(3) represents a particular (possibly hypothetical) disease. With respect to this particular disease, we also assume that the parameter values given in Table 1 represent some “nominal values”, which can be defined as the values that leads to a “best fit” of the model with the experimental data. We define “nominal model patient” to be a patient whose immunological dynamics conform to this mathematical model having these nominal parameter values (in addition, the model of a nominal model patient does not have any unwanted inputs). On the other hand, we define “realistic model patients” as those whose immunological dynamics conform to this model, but may have different values of the parameters. In addition, it may contain modelling error terms (e.g. unwanted inputs). For the rest of the paper, however, we use the terms “nominal patient” and “realistic patient” for better readability, assuming that they refer to “nominal model patient” and “realistic model patient” respectively. Even though the actual values for a real-

istic patient are expected to be around these nominal values, a major difficulty of a control system designer arises from the fact that these actual values are usually unknown, (it is practically impossible to come up with the exact values of these parameters for each individual patient before initiating the treatment). Hence for a control design strategy to be successful, it must be sufficiently robust to this issue of parameter uncertainty. In addition, it should also be robust to any “unknown” external disturbance (like an unknown rate of continued introduction of pathogens in the process under study).

2.2.

Objective of control design (medication strategy)

The objective of this paper is to present a system theoretic technique to come up with an adaptive drug delivery strategy (time history of the required drug delivery) by taking the help of neural networks so as to kill the invading microbes. At the same time the goal is also to make sure that the affected organ is not fatally damaged at any point of time during the treatment. In system theory language, the control synthesis strategy for designing an appropriate Ud (t) must make sure that [x1d , x4d ] → 0 as t → ∞ (i.e. as the treatment progresses). This objective should be met without any undesirable behavior in all the states, including the ones for which no specific goal is enforced. Note that meeting the objective as t → ∞ is written only to have a mathematically meaningful formulation. In practice (which will be clear in Section 4), the objective will be met much before t → ∞ (i.e. within a meaningful finite time). However, this mathematical formulation allows us to apply the drug administration scheme for an arbitrarily long amount of time, and hence, the duration of the drug administration need not be decided apriori. As mentioned before, achieving the above goal implies achieving only part of the overall objective. In general, actual system dynamics (which represent realistic patients) will have some parameter uncertainty (i.e. it will have some off-nominal parameter values) and possibly external disturbances as well. We denote the actual dynamics having these effects by ˙ = f (X) + B U +  X

(4)

where X = [x1 x2 x3 x4 ]T , U = [u1 u2 ]T are the state and control vectors respectively. The matrix B is as defined in (2) and the function f (X) is defined as

⎡

⎤

(a11 − a12 x1 )x3 a21 (x4 )a22 x1 x3 − a23 (x2 − x2∗ ) ⎥ ⎢ f (X)  ⎣ ⎦ a31 x2 − (a32 + a33 x1 )x3 a41 x1 − a42 x4

(5)

In an analogous manner, we assume that the parameter a21 is a nonlinear function of x4 given by

 a21 (x4 ) =

cos(x4 ), 0 ≤ x4 ≤ 1/2 0, x4 > 1/2

(6)

The numerical value of the parameter aij need not be same as its corresponding nominal value aijd . More impor-

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 4 ( 2 0 0 9 ) 207–222

tantly, in general the actual values of these parameters are “unknown”, even though their values are supposed to be around their corresponding nominal values.  represents the vector of unwanted inputs to the system, which is assumed to be bounded, smooth and slowly varying (so that it can be captured by a neural network trained online). While our mathematical formulation can tackle any such unwanted input, in our numerical studies we have confined ourselves to only a constant unknown signal that represents a constant rate of continuing introduction of pathogens. This was motivated from the analysis and numerical results presented in [35] and is sufficient to demonstrate the capability of the method presented. The goal of the nonlinear adaptive control design (treatment strategy) can now be defined as to come up with an appropriate U(t) to ensure that the actual states [x1 , x4 ] → 0 as t → ∞ (i.e. as the treatment progresses). This goal should be achieved despite the presence of parameter uncertainties and unwanted signal inputs. We meet this objective in two steps: (i) first we design a nominal controller Ud (t) that makes sure that [x1d , x4d ] → 0 as t → ∞; (ii) next, we design an adaptive controller U(t) (the adaptation being carried out online), which makes sure that [x1 , x4 ] → [x1d , x4d ] as t → ∞. This strategy makes the control design adaptive to parameter uncertainties as well as to unwanted inputs in the sense that we are guaranteed to have [x1 , x4 ] → 0 as t → ∞. Hence, the medication scheme is more effective in real-life. The mathematical details of the two steps are discussed in detail in Section 3. Once again, note that meeting the objective as t → ∞ is written only to have a mathematically meaningful formulation. In practice, the objective is met within a meaningful finite time.

3. Drug delivery scheme (control design): mathematical details In this section, we discuss the mathematical details of the drug delivery scheme, i.e. we outline the mathematical concepts of the control design procedure. However, note that the techniques used are not limited to the problem in discussion alone (rather they are applicable to a wide class of similar nonlinear problems). Hence we first discuss the relevant mathematics in a generic manner and then discuss the specific equations related to the particular drug delivery design problem.

3.1.

Nominal control design

As pointed out before, we first design a nominal controller which will meet the goals for the nominal model. Note that the neuro-adaptive technique proposed and used here (see Section 3.2) does not depend on any particular technique used to design the nominal controller. In this work we have used we have used ‘dynamic inversion’ technique for this purpose dynamic inversion mainly because of its elegance and simplicity. Even though dynamic inversion is a fairly established technique by now [4,14–18,32], we still outline the basic steps that are relevant to our problem for completeness of this paper.

3.1.1.

211

Generic theory (dynamic inversion)

We focus on a class of nonlinear systems which are affine in control and are represented by ˙ d = fd (Xd ) + Gd (Xd )Ud X

(7)

Yd = hd (Xd )

(8)

where Xd ∈ n , Ud ∈ m , Yd ∈ p are the state, control and performance output vectors of the nominal system respectively. We assume that m = p, i.e. we concentrate on ‘square systems’ for which the number of performance outputs are same as the number of control inputs. We also assume that the system is pointwise controllable [7]. The objective is to design a controller Ud so that Yd → Yd∗ as t → ∞, where Yd∗ (t) is the commanded signal for Yd to track. Yd∗ is assumed to be bounded, smooth and slowly varying. First using the chain rule of derivative, the expression for Y˙ d can be derived from (8) as Y˙ d = fYd (Xd ) + GYd (Xd )Ud

(9)

where fYd  [(∂hd /∂Xd )]fd (Xd ) and GYd  [(∂hd /∂Xd )]Gd (Xd ). Here we assume that GYd is not singular ∀t. Next, defining Ed  (Yd − Yd∗ ), the aim is to design Ud such that the following stable linear error dynamics is satisfied E˙ d + Kd Ed = 0

(10)

where Kd is chosen to be a positive-definite gain matrix (a relatively easier way is to choose it as a diagonal matrix with positive elements in the diagonal). For better physical interpretation, one can choose Kd = diag(1/1 , . . . , 1/m ) where i (i = 1, . . . , m) represents ‘time constant’ of the ith error channel. Using the definition of Ed , substituting the expression for Y˙ d from (9) in (10) and carrying out the necessary algebra we get the solution for the control variable as Ud = −[GYd (Xd )]

−1

{fYd (Xd ) + Kd (Yd − Yd∗ ) − Y˙ d∗ }

(11)

Before proceeding further, we wish to mention a few salient points with respect to this technique. First, note that it leads to a closed-form solution for the controller, and hence it can be implemented online without any computational difficulties. Moreover as long as (11) is satisfied, Ed → 0 ‘asymptotically’ (rather exponentially). In other words, asymptotic tracking is achieved. Despite these advantages, there are a few important issues with respect to this control design technique as well, which are outlined below: (i) The assumption m = p (same number of controls as outputs) need not hold good. If m < p, no perfect tracking is possible [20]. On the other hand, if m > p, additional objectives are usually introduced in the problem objective to obtain a meaningful solution for the controller [27]. However in our problem m = p, and hence, this issue does not arise. (ii) Another important issue is the question of ‘internal dynamics’. This concept is linked to ‘relative degree’  of the problem, which is defined as the total number

212

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 4 ( 2 0 0 9 ) 207–222

of derivatives taken in (9)(a first-order error dynamics in each output error channel is not always possible). Unless  = n, one explicitly needs to answer the question of internal stability [37]. Unless internal dynamics is stable, the synthesized controller is of no practical use. This issue does arise in our problem since in our case  = 2 and n = 4. Hence, we have explicitly analyzed this issue and have proved that the internal dynamics remains stable. (iii) For successful implementation of any control design, it must be demonstrated that it is sufficiently robust to the issues of parameter uncertainties in the model. Dynamic inversion approach, however, is quite sensitive to this issue in general [4,15,18]. This is also an issue of concern in the drug delivery design problem discussed in this paper, for which we have proposed to use the neuroadaptive method (discussed in Section 3.2).

3.1.2.

of parameter inaccuracies as long as signs of a23 and a32 (corresponding off-nominal values of a23d and a32d ) are same as those of a23d and a32d respectively (i.e. they are both positive).

3.2.

As pointed out before, for the successful implementation of any control strategy the issues of parameter uncertainty, modeling inaccuracies and unwanted inputs must be addressed. Here we present a model-following adaptive control design using neural networks, which is capable of addressing all of the above issues. Note that even though we have assumed no modeling inaccuracies explicitly, the difference of parameter values from their nominal values essentially generate unknown algebraic terms in the model. These unknown functions are captured by the neural networks (which are trained online) and subsequently the control solution is updated so that the state trajectories of the actual system follow the corresponding state trajectories of the closed loop nominal system. Once again since the mathematical development is generic enough to be applied to a wide class of nonlinear problems, we organize the mathematical details by outlining the generic theory first, which is followed by the specific equations with respect to the drug delivery problem under study.

Problem specific equations

With regard to the problem of treatment of infectious diseases, from (9) we select the performance output as Yd = [x1d x4d ]T and write the output dynamics as

Y˙ d =

(a11d − a12d x3 )x1 a41d x1 − a42d x4



+

b1 0

0 b4



u1d u2d

 (12)

3.2.1.

Note that in our case m = p = 2, fYd (Xd ) = T and GYd (Xd ) = [(a11d − a12d x3d )x1d (a41d x1d − a42d x4d )] diag(b1 , b4 ). Note also that GYd (Xd ) is a diagonal matrix with non-zero constants in its diagonal (therefore its inverse always exists). Since the goal is to make sure Yd = [x1d x4d ]T → 0, Yd∗ = 0 for our problem. Next, choosing Kd = diag(1/1d 1/2d ), substituting the required expressions in (11), we get the solution for the control variable (nominal controller) Ud as



u1d u2d





0 b4

b1 =− 0

−1

(a11d − a12d x3d )x1d a41d x1d − a42d x4d



+

1/1d 0

0 1/2d

Note that the relative degree for our problem  = 2, whereas the number of states n = 4. Hence, it is essential to prove the stability of internal dynamics. However, it is wellknown that it is sufficient to show that the “zero dynamics” is stable [37]. In our problem, the zero dynamics is essentially the homogeneous dynamics of the state variables x2d and x3d after Yd = [x1d x4d ]T → 0. To prove this, first we define the deviation variables x2d  x2d − x2∗ and x3d  x3d − x3∗ . Next, when d

d

Yd = [x1d x4d ]T = 0 we carry out the necessary algebra from (1) and (2) and derive the following deviation dynamics:



x˙ 2d x˙ 3d





= [Ad ]

x2d x3d





,

[Ad ] 

−a23d a31d

0 −a32d



Clearly, (14) turns out to be a linear system. It is well known from linear system theory [13] that this dynamics is asymptotically stable (rather exponentially stable) if and only if the eigenvalues of matrix Ad are in the left half of the complex plane. However, since Ad is a triangular matrix, its eigenvalues are its diagonal elements, which are −a23d = −1 and −a32d = −1.5. Hence we conclude that the internal dynamics is stable. Note that this stability is assured even in the presence

Neuro-adaptive design

Generic theory

Our focus in this paper is a class of nonlinear systems that are affine in control, for which the actual system dynamics is written as



˙ = f (X) + G(X)U +  X

(14)

Y = h(X)

(15)

x1d x4d

 (13)

where X ∈ n , U ∈ m , Y ∈ p (m = p) are the state, control and performance output vectors of the actual system respectively and  represents the vector of unwanted input to the system and it is assumed to be bounded, smooth and slowly varying (so that it can be captured by a neural network being trained online). Using chain rule of derivative in (15), the expression for Y˙ can be written as Y˙ = fY (X) + GY (X)U + Y

(16)

where fY  [∂h/∂X]f (X), GY  [∂h/∂X]G(X) and Y  [∂h/∂X]. Note that the functions fY and GY need not have the same algebraic expressions as fYd and GYd . We assume that the dynamics in (16) can be rewritten as Y˙ = fYd (X) + GYd (X)U + d(X)

(17)

where d(X) is an unknown function that arises due to parameter uncertainties, modeling inaccuracies and unwanted inputs. Functions fYd (X) and GYd (X) are similar to fYd (Xd ) and GYd (Xd ) respectively (they have the same parameter values). The objective of the problem is to design a suitable adaptive

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 4 ( 2 0 0 9 ) 207–222

213

controller U(t) such that Y(t) of (16)(which is same as that of (17)) tracks Yd (t) of (9); i.e. Y → Yd as t → ∞. In other words, the output of the actual system will be forced to track the output of the nominal model. The philosophy of this adaptive control design is to capture this unknown function d(X) through a neural network ˆ approximation d(X) and then use the updated model to compute the actual control U. However, the approach adopted here to ensure Y → Yd as t → ∞ involves two steps: (i) Y → Ya (an auxiliary variable) and (ii) Ya → Yd . These steps are discussed below. First we define the dynamics of the auxiliary variable Ya as follows: ˆ Y˙ a = fYd (X) + GYd (X)U + d(X) + K (Y − Ya )

(18) Fig. 1 – Structure of ith neural network.

where K > 0 (a positive definite matrix) can be selected as a diagonal matrix with the ith element being ki > 1. In the implementation of the proposed technique, we will also need the initial conditions, which are assumed as Ya (0) = Yd (0) = Y(0). The reason for choosing an auxiliary system dynamics of the form in (18) is to facilitate meaningful bounds on the tracking errors and neural network weights, which will be clear later in this section. First we discuss the philosophy of ensuring Ya → Yd as ˆ t → ∞, assuming that a neural network approximation d(X) of the unknown function d(X) is available to us. To achieve this objective, we propose to design the controller U such that the following stable error dynamics is enforced. (Y˙ a − Y˙ d ) + KY (Ya − Yd ) = 0

(19)

where the gain matrix KY is chosen as a positive definite matrix. Substituting (9) and (18) in (19), and carrying out the algebra, we get the solution for the control variable as U = −[GYd (X)]

−1

ˇ

(20)

where

y˙ i = fYd (X) + [GYd (X)U]i + di (X)

(23)

y˙ ai = fYd (X) + [GYd (X)U]i + dˆ i (X) + ki eai

(24)

i

i

Subtracting (24) from (23) and using the definition of eai we can write e˙ ai = di (X) − dˆ i (X) − ki eai

ˆ ˇ  fYd (X) + d(X) + K (Y − Ya ) − fYd (Xd ) − GYd (Xd )Ud + KY (Ya − Yd )

of linear-in-the-weight neural networks (Fig. 1), where the weights of the network appear linearly [8]. Note that in this formulation, the combination of psub-networks can be interpreted to constitute a single neural network that represents d(X). The idea of having p neural networks for p independent channels is to facilitate simpler mathematical analysis. More important, it leads to faster training as none of the weights are linked to more than one output function (which reduces computational complexity). Next, we present a way of updating the weights of the neural network (i.e. training them). In this regard, we first define the error eai  (yi − yai ). Then from (17) and (18), the actual and auxiliary dynamics of the ith channel can be written as

(21)

(25)

However, from the universal function approximation property of neural networks [8] we know that there exists an ideal neural network with an optimum weight Wi (with respect to the basis function ˚i (X)) such that it approximates di (X) with accuracy εi . In other words, we can write

Next, we discuss how to capture the unknown function ˆ d(X) as a neural network approximation d(X) which is crucial for the control solution in (20). In the process, we will ensure that Y → Yd as t → ∞. For this purpose, we first write T d(X)  [d1 (X)· · ·dp (X)] where di (X), i = 1, 2, . . . , p is the ith component of d(X). Next, we propose to approximate each di (X) as dˆ i (X) in a separate sub-neural network and write

˜ T ˚i (X) + εi − k ea e˙ai = W i i i

ˆ T ˚i (X) dˆ i (X) = W i

ˆ i ) is the error between ideal weight and ˜ i  (Wi − W where W ˜˙ i = −W ˆ˙ i actual weight of the neural network. Note that W

(22)

ˆ i ]h×1 is the weight vector of the ith neural network where [W and [˚i (X)]h×1 is its basis function vector (which is chosen judiciously by the control designer depending on the problem). The essential idea of writing (22) is to follow the principle

di (X) = Wi T ˚i (X) + εi

(26)

On substituting (22) and (26) in (25) we get (27)

(since Wi is a constant). Next, we define a series of positive definite Lyapunov function candidates Li , i = 1, . . . , p as follows: Li =

1 1 ˜ TW ˜ i) (W (ea ea ) + i 2 i i i 2i

(28)

214

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 4 ( 2 0 0 9 ) 207–222

where i , i > 0, ∀i = 1, . . . , p. Taking the time derivative of both sides of (28), we obtain ˜ TW ˜˙ i = ea i [W ˜˙ i ˜ T ˚i (X)+εi −k ea ]+ W ˜ T  −1 W L˙ i = eai i e˙ ai +i−1 W i i i i i i i ˜ T [ea i ˚i (X) −  −1 W ˆ˙ i ] + ea i εi − k e2a i =W i i i i i i

Next, the auxiliary output dynamics can be written as





(29)

Note that our objective is to come up with a meaningful condition that will ensure L˙ i < 0 so that we can use Lyapunov stability theory [14,32]. However, the expression for L˙ i contains ˜ i , which is unknown. Hence, nothing can be said about the W sign of L˙ i . To address this difficulty, we force the term multi˜ i to zero and obtain the following weight update rule plying W (training algorithm) for the neural network.



+



x1 − x1a x4 − x4a

+



0 b4

b1 0

u1 u2



+

dˆ 1 (X) dˆ 2 (X)



(34)

Assuming KY = diag(1/1Y 1/2Y ), where 1Y and 2Y are the time constants for the two error channels, substituting (12) and (34) in (19) and carrying out the necessary algebra in (20) we get the adaptive controller as



u1



 =−

u2

ˆ˙ i = i ea i ˚i (X) W i

(a11d − a12d x3 )x1 a41d x1 − a42d x4

Y˙ a =

b1

0

0

b4

−1 (ˇ1 − ˇ2 )

(35)

(30) where

At this point, it can be seen that the weight update dynamics involves a set of simple differential equations, which need to be integrated online using some numerical integration technique (e.g. forth-order Runge–Kutta technique [6]), which can be done easily without any computational complexity. For initial condition of the weights, typically it is assumed that ˆ i = 0. This is compatible with the argument that if the model W is perfect (i.e. there is no error), then dˆ i (X) = 0∀t. Note that i can be interpreted as the learning rate for the ith network and its numerical value essentially dictates the rate of convergence of dˆ i (X) to di (X). With (30) in place, (29) reduces to L˙ i = eai i εi − ki e2ai i

(31)

From (31) it is clear that L˙ i < 0 if e2ai > εi eai /ki . This eventually leads to the following condition: L˙ i < 0

whenever

|eai | > |εi /ki |

 ˇ1 

 ˇ2 

=

−

0

0

b4

u1d

 

u2d

1/Y1

0

0

1/Y2

−



(a11d − a12d x3 )x1d a41d x1 − a42d x4d

x1 − x1a



+

dˆ 1 (X)



dˆ 2 (X) (36)

x4 − x4a x1a − x1d

 



x4a − x4d

Next, we discuss the process of neural network selection and its training. The first step in this regard is to choose an appropriate basis function vector, for which we carry out the following analysis. Noting that the actual parameters in (33) are a result of uncertainties in the corresponding nominal parameters (we assume that there are no missing algebraic terms in the model), we can write



(32)

((a11d + a11 ) − (a12d + a12 )x3 )x1



(a41d + a41 )x1 − (a42d + a42 )x4

 +

 =

b1

0

0

b4



u1

 

u2

(a11d − a12d x3 )x1 a41d x1 − a42d x4

 +

+

ı1 0

  +



0

b1 0

a11 x1 − a12 x1 x3 + ı1



u1



u2

b4



(37)

a41 x1 − a42 x4

Problem specific equations

The output of the actual system is Y  [x1 x4 ]T . As suggested in [35], the external unwanted unknown input is the rate of continuing introduction of pathogen. So we can write Y = [ı1 0]T . Then from (16) and (17) we can write





b1

−

Y˙ =

Y˙ =

 

a41d x1 − a42d x4



Using Lyapunov stability theory [14], we conclude that the ˜ i will be pulled towards the origin trajectory of eai and W ˜ (eai = 0, Wi =0) as soon as |eai | > |εi |/ki . Hence, it is clear that the dynamics of error (27) and the weight update dynamics (30) are both “practically stable” [18]. Note that by appropriate selection of the size of the neural network and the basis functions, the value of i can be made small. Further more, as long as ki is chosen to be sufficiently high, the error bound will be very small for all practical purposes.

3.2.2.

(a11d − a12d x3 )x1

(a11 − a12 x3 )x1 a41 x1 − a42 x4



(a11d − a12d x3 )x1 a41d x1 − a42d x4

+

b1 0



+

b1 0

0 b4



0 b4

u1 u2





u1 u2

+

ı1 0



+



d1 (X) d2 (X)

 (33)

Comparing (37) with (33) it is clear that



d1 (X) d2 (X)



=

a11 x1 − a12 x1 x3 + ı1 a41 x1 − a42 x4

 (38)

Since as in [35] in our numerical simulations we assumed ı1 to be constant, we selected the basis functions for our neural networks as ˚1 (X) = [1 x1 x1 x3 ]T and ˚2 (X) = [x1 x4 ]T . It is important to note, however, that even though only a few parameters (namely a11 , a12 , a41 and a42 ) appear in (38),

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 4 ( 2 0 0 9 ) 207–222

in our simulation studies we have assumed uncertainties in all parameters and the neuro-adaptive technique was found to be working successfully. This is perhaps because of the generalized nature of the nueral network function approximation.

3.2.3.

Implementation of controller (medication scheme)

Even though the mathematical details presented for the completeness the control design process scheme seems to be a bit involved, we wish to point out that implementation of the algorithm developed is fairly straightforward. Since we have assumed a state feedback controller, the first step would be to sense the states of the system (namely the concentration of pathogens, plasma cells, antibodies and relative characteristic of the organ’s health). Since sensing all the states may not be possible, an observer (or better, a filter) design may be necessary for estimating the states, which is a topic of research by itself and is beyond the scope of this study. A first step in this direction can possibly be to use an extended Kalman filter with the dynamic inverse controller, as done in [33]. However, in this study we assume that the states are available for control computation. We assume also that the control update interval coincides with the step size chosen t and write the discrete time steps (at which controller is updated) as tk , k = 1, 2, . . .. With the above assumptions, the following steps would be required for implementing the proposed controller: (1) At initiation of the treatment (i.e. at t = t1 = 0), obtain the initial value of the state vector X1 = X(0) of the patient through diagnosis (we assume here that such a diagnosis process is available). Initialize Xd1 = Xa1 = X1 . Also initialˆ1 =W ˆ 2 = 0. ize the weights of the neural network as W 1 1 Next, at every instant of time tk , k = 1, 2 . . ., carry out the following steps: (2) Obtain the desired nominal controller Udk from (13). ˆ (3) Compute d(X) from (22) and obtain the controller Uk from

(4)

(5) (6) (7)

(35). This is the actual drug delivery that needs to be given to the patient at tk . Propagate the nominal and auxiliary states to tk1 (Xdk+1 and Xak+1 ) by integrating (12) and (34) by one time step t. This can be done by implementing an appropriate numerical integration scheme (we have used fourth-order Runge–Kutta method [6] for this purpose). In parallel, propagate the weights of the neural networks ˆ 1 ,W ˆ2 W by integrating (30). k+1 k+1 Obtain the value of the actual state Xk+1 (at tk+1 ) by diagnosis. If the patient is cured (which can be inferred from Xk+1 ), stop medication. Else go to (ii).

Note that even though for our simulation purpose we have assumed some numerical values of the actual parameters and integrated the corresponding state equation to obtain the actual state value Xk+1 at tk+1 , in real-life implementation this step is not necessary. Instead, the value of Xk+1 is supposed to be obtained directly from the patient (through continuous diagnosis). Hence, the knowledge of actual parameters and the value of external noise are not necessary. Here we reiterate that we have proposed this implementation with the

215

assumption that proper diagnosis procedures along with the necessary equipments are available to get the information about all states accurately without any substantial time-delay in the diagnosis process.

4.

Numerical results

4.1.

Treatment strategy for nominal patients

First, we present the results for the nominal system. For this purpose, we selected the initial condition for our numerical experiments based on the analysis carried out in [34] which can be summarized as follows. For a given set of parameter values, the homogenous system responds in four possible ways depending on the initial condition. These are represented by the sub-clinical, clinical, chronic and lethal cases. The external treatment becomes necessary for the chronic and lethal cases only as in the sub-clinical and clinical cases the body’s immune system is capable to cure the disease. As in [34], we have selected the lethal case for our numerical experiments for which the initial condition is given by T X(0) = [3 2 (a31d /a32d )x2∗ 0] . It is important to note that the technique presented in this study does not depend on the selection of initial condition (it works for any initial condition). However, we have presented the results with the selection of this initial condition for better physical meaning of the results. An interested reader can see the reference for more justification of this selection. The control design parameters were chosen as 1d = 0.5 and 2d = 1.5. Forth-order Runge–Kutta method with constant step size t = 1/60 was used for numerical integration of differential equations. Assuming that the time units of the generic model used is in hours, the selection of t = 1/60 implies that the control is updated every minute (biological systems, in general, are slow reacting system). It can be seen from Fig. 2 that the goal of [y1d y2d ] = [x1d x4d ] → 0 as t → ∞ is met without any problem. The concentration of pathogen (x1d ) decreases asymptotically to zero, which is expected from a dynamic inversion controller. Also, the organ remains healthy throughout the duration of treatment. The stability of internal dynamics is clearly observed in Fig. 2 as the concentrations of the plasma cells and antibodies are clearly bounded. Moreover, as expected, they tend towards their respective steady state values with time. The associated control histories are shown in Fig. 3. It can be observed that their magnitudes decay smoothly towards zero as well. Note that the medication can be stopped after t = 2units, since for all practical purpose the values of controls are zero thereafter.

4.2.

Treatment strategy for realistic patients

As pointed out before, even though an effective control therapy has been obtained for the nominal patient, it is seldom the case that realistic patients will have the same nominal parameters as used in the model. Rather, in most of the cases the parameter values will not match exactly as their corresponding nominal values. Before implementing the neuro-adaptive design, we wanted to carry out a robustness

216

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 4 ( 2 0 0 9 ) 207–222

Fig. 2 – State trajectories of the ideal patient.

simulation study (with respect to parameter uncertainties) and experiment with some off-nominal parameter values in the system. We perturbed all the parameter values within the function f (X) in the model (see Eq. (5)) and selected numerical values of the parameters randomly within 20% of their corresponding nominal values. In other words, we selected random values of parameters from the following ranges: a11 ∈ [0.8, 1.2], a12 ∈ [0.8, 1.2], a22 ∈ [2.4, 3.6], a23 ∈ [0.8, 1.2], a31 ∈ [0.8, 1.2], a32 ∈ [1.2, 1.8], a33 ∈ [0.4, 0.6], a41 ∈ [0.4, 0.6] and a42 ∈ [0.8, 1.2] and carried out the simulation studies. It was interesting to observe that for some cases the resulting drug delivery (controller) still works in the sense that it is able to cure the disease. This indicates that the synthesized controller based on nominal parameters has some inherent robustness. However, in some other cases it does not work, indicating that in such a case the patient will not be cured despite being on continuous medication. One such set of parameter values is given in Table 2. Note that the initial T condition was chosen as X(0) = [3 2(a31 /a32 )x2∗ 0] .

Table 2 – Actual parameter values. Parameter

Value

a11 a12 a21 a22 a23 a31

1.19 0.8116 See (6) 3.4821 0.8243 1.1783

Parameter a32 a33 a41 a42 b1 b4

Value 1.3524 0.5456 0.5019 0.9643 −1 −1

First we present the results when the nominal controller (medication) is applied to this realistic patient having offnominal parameter values (we assume that this patient is not under continuous pathogen attack). The results are plotted in Fig. 4. It can be seen that this treatment strategy is ineffective. The transient values of all the state variables (especially the organ health and concentration of pathogen) grow to unacceptably high values. Besides they also show unacceptably high oscillations. Most important of all, it can be

Fig. 3 – Drug delivery histories for the ideal patient.

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 4 ( 2 0 0 9 ) 207–222

217

Fig. 4 – State trajectories of an actual patient (with nominal drug delivery).

observed that the organ health indicator crosses unity, implying the death of the organ! Note that the graph plotted beyond the point where x4 crossing unity is for the sake of illustration alone and has no physical meaning. Furthermore, even if we assume a hypothetical case where the organ can somehow be recovered, the trajectories do not evolve towards their corresponding nominal steady state values of zero. This implies a chronic illness, where the patient is never completely cured. The corresponding treatment histories are shown in Fig. 5, which show bizarre behavior as well. It is obvious that the negative values of the pathogen killer (u1 ) can not be implemented. Most likely, it will not be possible to implement high values of the organ healer (u2 ) either. In any case, even if we assume that these drug delivery histories can be implemented, they are of no use as the strategy does not work. Next, we present the numerical results with the same parameter values (for which the nominal controller fails) by

applying the proposed neuro-adaptive controller. The control design parameters were chosen as 1Y = 2Y = 0.4, k1 = k2 = 1, 1 = 2 = 1, 1 = 1.2 and 2 = 10. Fig. 6 shows trajectories of the actual states on application of the modified control strategy. We observe that the neuro-adaptive controller was able to kill the invading microbes and preserve organ health, hence curing the hypothetical patient. This trend was observed in all of the numerous simulation studies we carried out. We have also compared the actual state trajectories with the corresponding state trajectories of the nominal model in Fig. 6. It is clear that the actual output [y1 y2 ]T = [x1 x4 ]T tracks the corresponding outputs [y1d y2d ]T = [x1d x4d ]T very well, which is quite compatible with the mathematical formulation. We also observe that even though no specific goal is enforced on the state variables x2 and x3 , these variables do not diverge away from the corresponding variables of the nominal model x2d and x3d . In Fig. 7, we have compared the control histories of the adaptive controller with the corresponding nominal (ideal)

Fig. 5 – Drug delivery histories for of an actual patient.

218

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 4 ( 2 0 0 9 ) 207–222

Fig. 6 – Comparison of nominal states using nominal drug delivery and actual states using adaptive drug delivery.

Fig. 7 – Comparison of nominal drug delivery for ideal patients and adaptive drug delivery for the actual patient.

Fig. 8 – Neural network approximation of unknown functions.

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 4 ( 2 0 0 9 ) 207–222

219

Fig. 9 – Histories of neural network weights.

Fig. 10 – Comparison of nominal states using nominal drug delivery and actual states using adaptive drug delivery.

controllers as applied to the nominal system. It is quite interesting to observe that the plots are quite close to each other. This observation, along with the fact that the output trajectories are also quite close to each other, implies that the neuro-adaptive formulation indeed enforces the output dynamics of the actual system to behave like that of the

ideal system and thereby makes the treatment strategy quite effective. Fig. 8 shows how the neural networks dynamically capture the unknown functions arising due to uncertainties in the parameters of the dynamic model. It can be seen that after sometime (i.e. after training), the networks capture the

Fig. 11 – Comparison of nominal drug delivery for ideal patients and adaptive drug delivery for actual patients.

220

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 4 ( 2 0 0 9 ) 207–222

Fig. 12 – Neural network approximation of unknown functions.

unknown functions quite well (with almost no error). Fig. 9 shows the training process of the weights of the neural networks. As expected, all the weights of the neural networks develop towards constant steady-state values. The trajectories of the weights in transition (i.e. while being trained) are also continuous and fairly smooth, which is an additional important benefit. Otherwise it will lead to fluctuations in drug delivery, which may either be uncomfortable for the patient or may even lead to undesirable side effects. So far we have not taken into account the effect of undesirable external disturbance inputs such as an unknown rate of introduction of pathogens. However, the neuro-adaptive technique presented here is capable of addressing the above issue directly, without any modification/augmentation of the basic technique. This is clearly a major advantage of the new approach. We demonstrate this claim in a set of representative numerical results. Even though the proposed technique should work for any unwanted input that is bounded and slowly varying, for demonstration purpose in our numerical simulations we have selected the unknown constant rate of introduction of pathogens as ı = 0.5 (same as [35]). Note that this is in addition to parameter uncertainties that we have accounted for (see Table 2). Once again we observed that the nominal controller fails to cure the patient. The plots corresponding to this failure case were observed to be quite similar to Figs. 4 and 5(and hence are not included here). However, when the neuro-adaptive controller is applied, the corresponding state and control histories are as shown in Figs. 10 and 11. From Fig. 10 it is clear that once again the modified treatment is effective in curing the disease by killing

the invading microbes and preserving organ health. It is also clear that the internal dynamics (i.e. concentrations of plasma cells and antibodies) is stable. In Fig. 11 the nominal and modified treatment histories are shown. As expected, they are also fairly close to each other. Note that the neuro-adaptive controller was successful in automatically finding the bias value of pathogen killer necessary to offset the effect of the unknown rate of pathogen attack. The neural network approximation of the unknown functions arising due to both parameter uncertainties (as in Table 2) as well as the external constant disturbance input is shown in Fig. 12. The history of the weights of the neural networks is shown in Fig. 13. Once again it can be seen that after training, the networks capture the unknown functions quite well (with almost no error) and all weights of the neural networks develop towards constant steady-state values. The trajectories of the weights in transition are also continuous and fairly smooth. In summary, all beneficial natures are preserved in this case as well. In addition, close observation of Fig. 13 as compared to Fig. 9 reveals the following. One weight of the first neural network (corresponding to y1 dynamics) in Fig. 9 develop towards zero in steady state, whereas it is a nonzero constant in Fig. 13. This weight history clearly demonstrates the capability of the proposed approach in successfully identifying a bias term. If there is no bias the corresponding weight automatically remain at zero. Before concluding this section, we would like to point out that even though only one set of results are presented in Figs. 7–11, similar successful results were obtained in many other cases as well. This clearly demonstrates that the pro-

Fig. 13 – Histories of neural network weights.

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 4 ( 2 0 0 9 ) 207–222

posed adaptive drug delivery design holds good promise as an effective way of treatment of the infectious diseases.

5.

Conclusions

Taking the help of advanced nonlinear system theory techniques of dynamic inversion and model-following neuroadaptive design, an automatic drug administration strategy is presented in this paper, which is applied to a generic model of the innate immune system. The method proposed in this paper adaptively adjusts the control (drug delivery) magnitudes online and is effective in killing the invading microbes and healing the damaged organ. It is adaptive with respect to parametric uncertainties and continued external disturbances. Simulation studies indicate that the adaptive controller proposed in this paper holds promise for its online implementation in the future. However it assumes that the required diagnosis process can be carried out at a sufficient faster rate so that all the states are available for control computation. The mathematical techniques presented are generic enough to be applied to other similar nonlinear control design problems as well. They are also computationally nonintensive and all required computations can be carried out online. Note that we have not used any filter design to estimate the states required for the control computation based on the output information (which can possibly be corrupted by sensor noise). However, such a study has many concerns of its own (like assuring filter stability, maximum noise filtering with minimum delay, studying the overall closed loop system performance, tuning the filter, etc.) and hence can be a topic for further research. A first step in this direction can possibly be to use an extended Kalman filter [33] with the proposed adaptive control design in this paper.

Conflict of interest statement None declared.

Acknowledgment The authors are very much thankful to Dr. R. Nayak, Formerly Professor at the Department of Microbiology and Cell Biology, Indian Institute of Science, Bangalore, India, for fruitful intellectual discussions on various aspects of the functioning of human immunology.

references

[1] A.L. Asachenkov, G.I. Marchuk, R.R. Mohler, S.M. Zuev, Immunology and disease control: a systems approach, IEEE Transactions on Biomedical Engineering 41 (10) (1994). [2] A. Asachenkov, G. Marchuk, R. Mohler, S. Zuev, Immunology and disease control: a systems approach, in: Disease Dynamics, Birkhauser, Boston, 1994.

221

[3] A.J. Calise, B-J. Yang, J.I. Craig, Augmenting adaptive approach to control of flexible systems, Journal of Guidance, Control and Dynamics 27 (3) (2004) 387–396. [4] D. Enns, D. Bugajski, R. Hendrick, G. Stein, Dynamic inversion: an evolving methodology for flight control design, International Journal of Control 59 (1) (1994) 71–91. [5] J.D. Farmer, N.H. Packard, A.S. Perelson, The Immune system, adaptation and machine learning, Physica 22D (1986) 187–204. [6] S.K. Gupta, Numerical Methods for Engineers, Wiley Eastern Ltd., 1995. [7] K.D. Hammet, C.D. Hall, D.B. Ridgely, Controllability issues in nonlinear statedependent Riccati equation control, Journal of Guidance, Control, and Dynamics 21 (5) (1998) 767–773. [8] M.H. Hassoun, Fundamentals of Artificial Neural Networks, MIT Press, Cambridge, MA, 1995. [9] N. Hovakimyan, F. Nardi, A.J. Calise, H. Lee, Adaptive output feedback control of a class of nonlinear systems using neural networks, International Journal of Control 74 (12) (2001) 1161–1169. [10] Z. Huang, S.N. Balakrishnan, Robust neurocontrollers for systems with model uncertainties: a helicopter application, AIAA Journal of Guidance, Control and Dynamics 28 (3) (2005) 516–523. [11] K.J. Hunt, R. Zbikowski, D. Sbarbaro, P.J. Gawthorp, Neural networks for control systems—a survey, Automatica 28 (6) (1992) 1083–1112. [12] C.A. Janeway, P. Travers, M. Walport, M. Shlomchik, Immunobiology, Garland, London, 2001. [13] T. Kailath, Linear Systems, Prentice Hall, 1980. [14] H.K. Khalil, Nonlinear Systems, 3rd ed., Prentice Hall Inc., NJ, USA, 2002. [15] B.S. Kim, A.J. Calise, Nonlinear flight control using neural networks, AIAA Journal of Guidance, Control and Dynamics 20 (1) (1997) 26–33. [16] D.E. Kirk, Optimal Control Theory: An Introduction, Prentice Hall, 1970. [17] S.H. Lane, R.F. Stengel, Flight control design using non-linear inverse dynamics, Automatica 24 (4) (1988) 471–483. [18] J. Leitner, A. Calise, J.V.R. Prasad, Analysis of adaptive neural networks for helicopter flight controls, AIAA Journal of Guidance, Control and Dynamics 20 (5) (1997) 972–979. [19] F.L. Lewis, A. Yesildirek, K. Liu, Multilayer neural net robot controller with guaranteed tracking performance, IEEE Transactions on Neural Networks 7 (2) (1996) 388–399. [20] Y. Li, N. Sundararajan, P. Saratchandran, Neuro-controller design for nonlinear fighter aircraft maneuver using fully tuned RBF networks, Automatica, 37 (8) (2001) 1293–1301. [21] P.M. Lydyard, A. Whelan, M.W. Fanger, Dynamic inversion: an evolving methodology for flight control design, in: Instant Notes in Immunology, Springer–Verlag, New York, 2000. [22] G.I. Marchuk, Some mathematical models in immunology, in: Proceedings of the IFIP Conference Optimization Techniques, Springer–Verlag, Heidelberg, 1978, pp. 41–62. [23] W.T. Miller, R. Sutton, P.J. Werbos, Neural Networks for Control, MIT Press, 1990. [24] R.R. Mohler, C. Bruni, A. Gandolfi, A systems approach to immunology, Proceedings of the IEEE 68 (1980) 964–990. [25] G.J.V. Nossal, Antibodies and Immunity, Basic Books, New York, 1969. [26] R. Padhi, N. Unnikrishnan, S.N. Balakrishnan, Model following robust neuro-adaptive control design for non-square, non-affine systems, in: ASME International Mechanical Engineering Congress and Exposition, 2004. [27] R. Padhi, M. Kothari, A neuro-adaptive augmented optimal dynamic inversion approach for effective treatment of

222

[28] [29] [30]

[31]

[32]

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 9 4 ( 2 0 0 9 ) 207–222

chronic myelogenous leukemia, in: Proceedings of the 2006 American Control Conference, June 14–16, 2006. A.S. Perelson, Immunology for physicists, Reviews of Modern Physics 69 (4) (1997) 1219–1267. I. Roitt, Essential Immunology, Oxford, 1974. R.M. Sanner, J.J.E. Slotine, Gaussian networks for direct adaptive control, IEEE Transactions on Neural Networks 3 (6) (March 1996) 837–863. S. Seshagiri, H.K. Khalil, Output feedback control of nonlinear systems using RBF neural networks, IEEE Transactions on Neural Networks 11 (1) (2000) 69– 79. J.J.E. Slotine, W. Li, Dynamic inversion: an evolving methodology for flight control design, in: Applied Nonlinear Control, Prentice Hall, 1991.

[33] R.F. Stengel, S.S. Mulgund, Optimal nonlinear estimation for aircraft flight control in wind shear, Automatica 32 (1) (1996) 3–13. [34] R.F. Stengel, R. Ghigliazza, N. Kulkarni, O. Laplace, Optimal control of innate immune response, Optimal Control Applications and Methods 23 (2002) 91–104. [35] R.F. Stengel, R.M. Ghigliazza, N.V. Kulkarni, Optimal enhancement of immune response, Bioinformatics 18 (9) (2002) 1227–1235. [36] R.F. Stengel, R. Ghigliazza, Stochastic optimal therapy for enhanced immune response, Mathematical Biosciences 191 (2004) 123–142. [37] E.M. Wallner, K.H. Well, Attitude control of a reentry vehicle with internal dynamics, AIAA Journal of Guidance, Control and Dynamics 26 (6) (2003) 846–854.