Chapter 1
Introduction to Partial Differentiation Equation Analysis: Chemotaxis

1.1 Introduction

This chapter serves as an introduction to the analysis of biomedical science and engineering (BMSE) systems based on partial differential equations (PDEs) programmed in R. The general format of this chapter and the chapters that follow consists of the following steps:

• Presentation of a PDE model as a system of PDEs, possibly with the inclusion of some additional ordinary differential equations (ODEs).
• Review of algorithms for the numerical solution of the PDE model.
• Discussion of a set of R routines that implement the numerical algorithms as applied to the model.
• Review of the computed output.
• Conclusions concerning the model, computer implementation, output, and possible extensions of the analysis.

This format is introductory and application oriented with a minimum of mathematical formality. The intention is to help the reader start with PDE analysis of BMSE systems without becoming deeply involved in the details of PDE numerical methods and their computer implementation (e.g., coding). Also, the presentation is self-contained so that the reader will not have to go to other sources such as a software download to find the routines that are discussed and used in a particular application. Our final objective then is for the detailed discussions of the various applications to facilitate a start in the PDE analysis of BMSE systems.

In this chapter, we consider the following topics.

• A brief introduction to PDEs.
• Application of PDE analysis to chemotaxis.
• Algorithms for the numerical solution of a simultaneous 2-PDE nonlinear chemotaxis model.
• Computer routines for implementation of the numerical algorithms.
• Traveling wave features of the 2-PDE chemotaxis model numerical solution.

1.2 Linear Diffusion Model

Inanimate systems have the general feature wherein chemical species move from regions of high concentration to regions of low concentration by mechanisms that are often modeled as diffusion, that is, according to Fick's first and second laws. In 1D, this diffusion is described (according to Fick's second law) as

where

In accordance with the usual convention for PDE notation, the dependent variable will subsequently be denoted as rather than . Thus, eq. (1.1a) will be

We can note the following features of eq. (1.1a).

• Eq. (1.1a) is a PDE because it has two independent variables, and . A differential equation with only one independent variable is termed an ordinary differential equation (ODE). Note also that is used to denote a partial derivative.
• The solution of eq. (1.1a) is the dependent variable as a function of the independent variables and , that is, in numerical form (rather than analytical form).
• Eq. (1.1a) is linear for constant D because the dependent variable and its partial derivatives are to the first degree (not to be confused with order because eq. (1.1a) is first order in because of the first-order derivative in and second order in because of the second-order derivative in ). Classifying eq. (1.1a) as linear presupposes that the diffusivity is not a function of . Eq. (1.1a) is nonlinear if because of the product .
• The diffusivity is inside the first (left most) differentiation to handle the case when is a function of and/or . If is a constant, it can be moved outside the first differentiation.

Eq. (1.1a) models ordinary diffusion because of, for example, random motion of molecules. A distinguishing feature of this type of diffusion is net movement in the direction of decreasing concentration as reflected in Fick's first law

is a component of the diffusion flux vector (with additional components and in Cartesian coordinates) and is therefore denoted with boldface. The minus sign signifies diffusion in the direction of decreasing concentration (as for ).

Since eq. (1.1a) is first order in and second order in , it requires one initial condition (IC) and two boundary conditions (BCs). For example,

where is a constant (length) to be specified and are functions to be specified. Since BCs (1.1d,e) specify the dependent variable at two particular (boundary) values of , that is, , they are termed Dirichlet BCs. The derivatives in can be specified as BCs, for example,

Eqs. (1.1f) and (1.1g) are termed Neumann BCs. Also, the dependent variable and its derivative can be specified at a boundary, for example,

Eqs. (1.1h) and (1.1i) are termed third-type, Robin, or natural BCs. All of these various forms of BCs (eqs. (1.1d)–(1.1i)) are useful in applications.

Finally, note that is defined over an open-ended interval or domain, , and is termed an initial value variable (typically time in an application). is defined between two different (boundary) values in , denoted here as or more generally (typically physical boundaries in an application). However, the interval in can be semi-infinite, for example, or fully infinite, .

1.3 Nonlinear Chemotaxis Model

Eqs. (1.1a) and (1.1b) can be extended to a nonlinear form of Fick's first and second laws. Also, we can consider more than one dependent variable, and we now consider two dependent variables and in place of just . The 2-PDE model for chemotaxis is ([2], p 68)

Here, we have employed subscript notation for some of the partial derivatives. For example,

Note that the PDE variables can have two subscripts. The first is a number denoting a particular dependent variable. The second is a letter denoting a partial derivative with respect to a particular independent variable. For example, denotes the first dependent variable differentiated with respect to . Also, the second (letter) subscript can be repeated to denote a higher order derivative. For example, subscript notation can be used in eq. (1.2b),

This compact subscript notation for partial derivatives can be useful in conveying a correspondence between the mathematics and the associated computer coding. This will be illustrated in the subsequent programming of eqs. (1.2).

The variables and parameters in eqs. (1.2) are

Eqs. (1.2) defines the volume concentration of a microorganism (such as bacteria), , when responding to an attractant (such as a nutrient or food supply), . Eq. (1.2a) reflects the rate of consumption of , that is, , due to ; the rate constant is taken as so that the minus is required for consumption ( with ).

Eq. (1.2b) is an extension of eq. (1.1b) and implies a diffusion flux in the -direction.

(the subscript in has been dropped). Eq. (1.2c) can be considered an extension of eq. (1.1b). It is nonlinear because of the term . We can note the following details about eq. (1.2c).

• The first RHS term, , is just Fick's first law, eq. (1.1b). In other words, the flux of eq. (1.2c), , is composed partly of the usual flux in the direction of decreasing gradient (, which tends to make because of the minus sign in eq. (1.2c)).
• The second RHS term is opposite in sign to the first and, therefore, gives the opposite effect for the flux . Note that the gradient in this term is , not . Thus, this term causes the bacteria flux to increase with increasing attractant concentration . Also, the ratio is a factor in determining the flux . This ratio causes the rate of transfer (flux) of the bacteria, , to increase with increasing bacteria concentration, , and also to increase with a decrease in the attractant concentration, ; for the latter, the bacteria apparently move faster when facing decreasing availability of the attractant (e.g., nutrient or food). Clearly, this nonlinear RHS term is a significant departure from the diffusion term of Fick's first law, eq. (1.1b). This is a unique feature of chemotaxis by which the bacteria seek higher concentrations of the attractant; this seems plausible if, for example, the attractant is a nutrient such as food. This effect is clearly a feature of an animate (living) system, such as bacteria, rather than an inanimate system.
• As the second RHS term has a rather unconventional form, we would expect that it will introduce unusual features in the solution when compared with the usual diffusion modeled by Fick's first and second laws. These features will be considered when the routines for the solution of eqs. (1.2a) and (1.2b) are discussed subsequently.
• As the nonlinear diffusion term ...