# Delay partial differential equations

Post-publication activity

Curator: Barbara Zubik-Kowal

A delay partial differential equation (DPDE) is an equation which involves

• at least two independent variables,
• an unknown function of the independent variables,
• the behavior of the unknown function at some prior value(s) of the independent variable(s),
• partial derivative(s) of the unknown function with respect to the independent variable(s).

Therefore, a delay partial differential equation differs from a partial differential equation in that it depends not only on the solution at a present stage but also on the solution at some past stage(s). If, additionally, the equation depends on the derivative(s) of the solution at some past stage(s), then it is a neutral delay partial differential equation.

Delay partial differential equations are also called partial functional differential equations as their unknown solutions are used in these equations as functional arguments.

## Applications of Delay Partial Differential Equations

Delay partial differential equations arise from various applications, like biology, medicine, control theory, climate models, and many others (see e.g. Wu (1996) and the references therein). Their independent variables are time $$t$$ and one or more dimensional variable $$x\ ,$$ which usually represents position in space but may also represent relative DNA content, size of cells, or their maturation level, or other values. The solutions (dependent variables) of delay partial differential equations may represent temperature, voltage, or concentrations or densities of various particles, for example cells, bacteria, chemicals, animals and so on.

Delay partial differential equations depend on their solutions at retarded variables or/and they depend on smooth averages of the solutions over some retarded intervals. Examples of such cases are given in the sections below with

• time delays $$\rightarrow$$ (1),
• time delay and shrinking in $$x$$ $$\rightarrow$$(2),
• integral over past time interval $$\rightarrow$$ (3),
• time delay, stretching in $$x\ ,$$ and integral over $$x$$-domain $$\rightarrow$$ (4)-(5),
• time delay in the derivatives of the solution $$\rightarrow$$ (6).

## Time Delays

Automatically Controlled Furnace

An example of delay partial differential equations is $\tag{1} \frac{\partial u(t,x)}{\partial t} = k \frac{\partial^2 u(t,x)}{\partial x^2} + v\Bigg(g\Big(u(t-\tau_1,x)\Big)\Bigg) \frac{\partial u(t,x)}{\partial x} + c\Bigg( f\Big( u(t-\tau_2,x) \Big) - u(t,x) \Bigg).$

Equation (1) is a parabolic delay partial differential equation with two delays $$\tau_1>0$$ and $$\tau_2>0\ .$$ It has been proposed by Wang (1963) to consider an automatically controlled furnace, see Fig.1 in Wang (1963). From this reference, the furnace is fed by the material strip that has to be heat-treated with a controlled temperature. The furnace temperature is varied by means of a heater actuated by a heater controller. The control objective is to maintain a desired spatial temperature distribution in the incoming material, which is fed into the furnace by rollers, the speed of which is regulated by a speed controller. This may be accomplished by placing temperature transducers along the material strip. The transducers provide information for a computer, which generates the appropriate control signals for the heater and feed-roller speed controllers.

In (1), $$x$$ means a position in the furnace, $$t$$ time, $$u(t,x)$$ a spatial average of the temperature distribution, $$v$$ is the instantaneous material strip velocity depending on the time-delayed temperature distribution $$u(t-\tau_1,x)\ ,$$ and $$f$$ represents a distributed temperature source function depending on $$u(t-\tau_2,x)\ .$$ The control signals may be delayed in time due to the possible presence of time delays in actuations, information transmission and processing. Therefore, the behavior of the overall control system depends not only on the present values $u(t,x),\quad \frac{\partial u(t,x)}{\partial x},\quad \frac{\partial^2 u(t,x)}{\partial x^2}$ but also on the time-delayed values $u(t-\tau_1,x), \quad u(t-\tau_2,x).$

Because of the delays $$\tau_1$$ and $$\tau_2$$ equation (1) has to be supplemented by a different initial condition (see Figure 1) than the initial conditions for non-delay partial differential equations.

## Time Delay and Shrinking in x

Proliferating and Maturing Cellular Populations

Another example of delay partial differential equations is $\tag{2} \frac{\partial u(t,x)}{\partial t} + r\frac{\partial}{\partial x}\Bigg( xu(t,x)\Bigg) = -\gamma u(t,x) + \lambda u\Big(t-\tau ,\alpha x\Big) \Bigg(1-u\Big(t-\tau ,\alpha x\Big)\Bigg).$

Equation (2) is a hyperbolic delay partial differential equation with delay $$\tau >0$$ and shrinked argument $$\alpha x$$ (here $$0<\alpha <1$$). This equation has been derived by Rey & Mackey (1993) and describes cell population dynamics, in which there is simultaneous proliferation and maturation. The unknown function $$u(t,x)$$ is a total number of proliferating cells at time $$t$$ and maturation level $$x\ .$$ Here, each cell is characterized by the maturation variable $$x\ ,$$ the range of which can be taken, without loss of generality, from the interval $$[0,1]\ .$$ The last term in (2) is a delay term as it depends on the population of cells at an earlier time $$t-\tau$$ and a maturation level $$\alpha x\ .$$

For this model, it is assumed that the age $$\tau$$ at cytokinesis is identical between cells. It is also assumed that cells die at a rate $$\gamma >0$$ independent of age or maturation level, and mature with a velocity, which is proportional to the maturity $$x\ ,$$ that is $$V(x)=rx\ ,$$ where $$r>0$$ is a constant of proportionality. Under these assumptions, $$\alpha =e^{-r\tau }$$ and $$\lambda =e^{(\gamma +r)\tau }$$ (see Rey & Mackey (1993)). Figure 2 illustrates the initial and boundary conditions for (2).

## Integral over Past Time Interval

Single-species population models

Another example of delay partial differential equations is $\tag{3} \frac{\partial u(t,x)}{\partial t} = D\frac{\partial^2 u(t,x)}{\partial x^2} + r u(t,x) \Bigg(1-\int_{-\tau }^0u(t+s,x)d\eta (s)\Bigg).$

Equation (3) is a parabolic delay partial differential equation with an integral over past time intervals $$[t-\tau ,t]\ .$$ It has been proposed by Green & Stech (1981) for a class of single-species population models, where species can diffuse. From Green & Stech (1981), the unknown function $$u(t,x)$$ represents a size of a population at a time $$t$$ and position $$x\ .$$ The constants $$D$$ and $$r$$ are positive and determine rates of diffusion and growth, respectively. The initial and boundary conditions for (3) are illustrated in Figure 3.

Equation (3) generalizes the model by Hutchinson (1948), see also Kuang (1993), which can be obtained from (3) by replacing the integral by a single retarded value $$u(t-\tau )$$ and by putting $$D=0$$ (no diffusion) or by a single retarded value $$u(t-\tau,x)$$ for $$D\neq 0$$ (Hutchinson equation with diffusion), see Wu (1996). More complex models than (3) have been proposed by e.g. Gourley et al. (2004) and Gourley & Kuang (2005), see also references therein.

## Time Delay, Stretching in x, and Integral over x-Domain

Cancer Cell Division Cycle

An example of a system of delay partial differential equations is $\tag{4} \left\{ \begin{array}{l} \frac{\partial }{\partial t}G_1(t,x)=4bM(t,2x)-(k_1+\mu_{G_1}) G_1(t,x), \\ \frac{\partial }{\partial t}S(t,x)= D\frac{\partial^2}{\partial x^2}S(t,x)-\mu_SS(t,x)- g\frac{\partial }{\partial x}S(t,x)+k_1G_1(t,x)-I(t,x;T_S), \\ \frac{\partial }{\partial t}G_2(t,x)=I(t,x;T_S)- (k_2+\mu_{G_2})G_2(t,x), \\ \frac{\partial }{\partial t}M(t,x)=k_2G_2(t,x)-bM(t,x)- \mu_MM(t,x), \end{array} \right.$

where the delay term $$I(t,x;T_S)$$ is given by $\tag{5} I(t,x;T_S)= \left\{ \begin{array}{ll} \int_0^{\infty}k_1G_1(t-T_S,y) \gamma(T_S,x,y)dy, \quad & t\geq T_S,\\ 0, \quad & t<T_S. \end{array} \right.$

System (4)-(5) has been proposed by Basse et al. (2003). It describes four distinct phases, namely the $$G_1$$-phase, DNA synthesis or $$S$$-phase, $$G_2$$-phase, and mitosis or $$M$$-phase, which are considered as parts of the cancer cell division cycle. Fig.1 in Basse et al. (2003) expresses the accumulation of cells in each of the phases and the movement of cells between them.

In (4)-(5), $$t$$ is time (measured in hours) and $$x$$ is the dimensionless relative DNA content. The dependent variables, $$G_1(t,x)\ ,$$ $$S(t,x)\ ,$$ $$G_2(t,x)$$ and $$M(t,x)$$ represent the density of cells in the corresponding phases. The parameters $$\mu_{G_{1}}\ ,$$ $$\mu_S\ ,$$ $$\mu_{G_2}\ ,$$ and $$\mu_M$$ are the death rates in $$G_1\ ,$$ $$S\ ,$$ $$G_2\ ,$$ and $$M$$-phases, respectively. The parameters $$k_1$$ and $$k_2$$are the transition probabilities of cells from $$G_1$$to $$S$$-phase and from $$G_2$$ to $$M$$-phase, respectively; $$b$$ is the division rate; $$D$$ is the dispersion coefficient; $$g$$ is the average growth rate of DNA in the $$S$$-phase; and $$\gamma(\tau,x,y)$$ is a given weight function (see Basse et al. (2003)). Figure 4 illustrates the initial and boundary conditions for (4).

For the case of dispersion (DNA contents vary), $$D\neq 0$$ and the second equation in (4) is a parabolic delay partial differential equation. For the case of no dispersion (constant DNA contents), $$D=0$$ and it is a hyperbolic delay partial differential equation.

## Neutral Delay Partial Differential Equation

Coupled Oscillators

An example of neutral delay partial differential equations is $\tag{6} \frac{\partial }{\partial t} Qu_t = K \frac{\partial^2 }{\partial t^2} Qu_t +f(u_t).$

Equation (6) has been proposed by Hale (1994a & b) to consider a continuum of diffusively coupled oscillators, where the constant $$K>0$$ represents the diffusive interaction between neighboring oscillators. (6) is a continuous version of a system of ordinary neutral delay differential equations investigated in Wu & Xia (1996) (see also the introduction to ordinary neutral delay differential equations) for oscillators on a linear periodic lattice. Hale (1994a & b) obtained (6), from the system investigated in Wu & Xia (1993), by taking the limit as $$h\to 0\ ,$$ where $$h$$ is a spacing between the oscillators.

In (6), the value $$u(t,x)$$ represents the voltage at time $$t\geq 0$$ and position $$x\in S^1$$ in the unit circle $$S^1\ .$$ For (6), Hale (1994a & b) has defined $Qu_t = u(t,x)-qu(t-\tau ,x),$ with $$0\leq q<1$$ and $$\tau >0\ .$$ Here, the lower index in $$u_t$$ means the solution segment defined by $\tag{7} u_t(\theta ,x) = u(t+\theta ,x), \quad -\tau \leq \theta \leq 0.$

Definition (7) is similar to the definition that is used by Hale (1977) for functions of one variable $$t\ .$$

## Partial Functional Differential Equations

Cell Population

Although delay partial differential equations are also called partial functional differential equations, the class of functional equations is wider than the class of delay equations. An example of a partial functional differential equation, which is not a delay partial differential equation, is the following $\tag{8} \frac{\partial u(t,x)}{\partial t} + \frac{\partial}{\partial x}\Big( g(x)u(t,x)\Big) = -\Big( \mu(x) + b(x)\Big) u(t,x) +4b(2x)u(t,2x).$

Equation (8) has been formulated by Diekmann et al. (1984). It describes the growth of a size-structured cell population reproducing by fission into two identical daughters. The independent variables $$t$$ and $$x$$ denote time and size of cells, respectively. The dependent variable $$u$$ is an unknown density function. The functions $$\mu\ ,$$ $$b\ ,$$ and $$g$$ are the rates at which cells of size $$x$$ die, divide, and grow, respectively.

In (8), none of the independent variables is retarded. However, (8) does contain a transformed $$x$$-variable, namely $$2x$$ in $$u(t,2x)\ ,$$ and $$(t,x)$$ is not the only argument in this equation. Therefore, $$u$$ is used in (8) as a functional argument and (8) is called a partial functional differential equation.

Other examples of partial functional differential equations are the delay partial differential equations (1), (2), (3), (4)-(5), and (6), where the unknown solutions act as functional arguments because $$(t,x)$$ is not the only point at which their behavior is needed for the equations. More complex models than (8) have been proposed by Hoppensteadt (1977).

## Initial and Boundary Conditions

Delay partial differential equations are supplemented by different kinds of initial and boundary conditions. For example, (1) is supplemented by $\tag{9} \left\{ \begin{array}{l} u(t,x)=u_0(t,x),\quad -\max\{ \tau_1,\tau_2\} \leq t \leq 0,\quad 0< x < 1,\\ u(0,t)=0,\quad u(1,t)=0, \quad t\geq 0, \end{array} \right.$

where $$u_0$$ is a known initial function. Here, the width of the furnace, which is considered in Wang (1963), is normalized and equal to $$1\ .$$ The conditions (9) are illustrated in Figure 1 by the blue rectangle (initial set) and the red segments (bounderies). The initial function $$u_0$$ is defined on the blue rectangle. Therefore, the values $$u(t_0-\tau_1,x_0)=u_0(t_0-\tau_1,x_0)\ ,$$ $$u(t_0-\tau_2,x_0)=u_0(t_0-\tau_2,x_0)$$ are known from (9) and can be used for (1) to determine $$u(t_0,x_0)\ .$$ The value $$u(t_1,x_1)$$ is joined in (1) with $$u(t_1-\tau_1,x_1)=u_0(t_1-\tau_1,x_1),$$ known from (9), and with $$u(t_1-\tau_2,x_1)\ ,$$ which is not defined by (9) but can be determined like $$u(t_0,x_0)\ .$$ The value $$u(t_2,x_2)$$ is joined in (1) with $$u(t_2-\tau_1,x_2)\ ,$$ $$u(t_2-\tau_2,x_2)\ ,$$ which are not defined by (9), and have to be determined from (1).

Equation (2) is supplemented by the initial-boundary conditions $\tag{10} \left\{ \begin{array}{l} u(t,x)=\varphi(x),\quad -\tau \leq t \leq 0,\quad 0< x < 1,\\ u(0,t)=\psi(t),\quad t\geq 0. \end{array} \right.$

Rey & Mackey (1993) have considered different initial and boundary functions $$\varphi$$ and $$\psi$$ for (2) and (10). Conditions (10) are illustrated in Figure 2 by the blue rectangle, for the initial condition, and the red segment for the boundary condition. Here, the value $$u(t_0-\tau,\alpha x_0)=u_0(t_0-\tau,\alpha x_0)$$ is known from (10) and can be used in (2) to determine $$u(t_0,x_0)\ .$$ On the other hand, $$u(t_1-\tau,\alpha x_1)$$ is not given by (10) and both values, $$u(t_1,x_1)$$ and $$u(t_1-\tau,\alpha x_1)$$ have to be determined from (2).

The initial and boundary conditions for (3) are written in the form $\tag{11} \left\{ \begin{array}{l} u(t,x)=u_0(t,x),\quad -\tau \leq t \leq 0,\quad 0< x < \pi,\\ u(0,t)=0,\quad u(\pi ,t)=0, \quad t\geq 0 \end{array} \right.$

(see Green and Stech (1981)). In Figure 3, the initial condition is illustrated by the blue rectangle and the boundary conditions are illustrated by the red segments. For (3), the integrals $$\int_{-\tau}^{0} u(t_0+s,x_0) ds$$ and $$\int_{-\tau}^{0} u(t_1+s,x_1) ds$$ have to be determined differently. Here, $\tag{12} \int_{-\tau}^{0} u(t_0+s,x_1) ds = \int_{-\tau}^{-t_0} u_0(t_0+s,x_1) ds + \int_{-t_0}^{0} u(t_0+s,x_1) ds$

and, by (11), the first integral on the right hand side of (12) can be determined from the known initial function $$u_0\ .$$ For the last integral in (12), the values of $$u$$ are not known and have to be determined from (3). Similarly, the values of $$u$$ for $$\int_{-\tau}^{0} u(t_1+s,x_1) ds$$ have to be determined from (3).

In accordance with experimental evidence, (4)-(5) are supplemented by the initial conditions $\tag{13} \left\{ \begin{array}{ll} G_1(0,x)=\frac{a_0}{\sqrt{2\pi \theta_0^2}} \exp\Bigg(-\frac{(x-1)^{2}}{2\theta_0^{2}}\Bigg), & 0<x<\infty, \\ S(0,x)=0, \quad G_2(0,x)=0, \quad M(0,x)=0, & 0<x<\infty, \end{array} \right.$

and the boundary condition $\tag{14} D\frac{\partial S(t,0)}{\partial x}-g\,S(t,0)=0, \quad t>0.$

Figure 4 illustrates (13)-(14). Note that since $$I(t,x;T_S)=0\ ,$$ for $$t<T_S$$ (by (5)), no values of $$G_1$$ are needed for $$t<0$$ and the initial set is the positive part of the $$x$$-axis (in blue). Therefore, the value $$S(t_0,x_0)$$ can be determined from (4) with $$I(t_0,x_0;T_S)=0\ .$$ However, $$S(t_1,x_1)$$ has to be determined from (4) with $$I(t_1,x_1;T_S)=\int_{0}^{\infty }k_1G_1(t_1-T_S,y)\gamma (T_S,x,y)dy$$ (by (5)), which needs information about the behavior of $$G_1$$ on the green line $$t=t_1-T_S\ .$$ (14) is the same zero flux boundary condition that often supplements non-delay partial differential equations.

Initial and boundary conditions for (8) are written in the form $\tag{15} \left\{ \begin{array}{l} u(0,x)=u_0(x),\quad \frac{a}{2}\leq x \leq 1\\ u(t,\frac{a}{2})=0,\quad t>0. \end{array} \right.$

They have been formulated by Diekmann et al. (1984) and are illustrated in Figure 5. Equation (8) is imposed for $$\frac{a}{2}\leq x \leq 1\ ,$$ $$t>0\ .$$ Here, $$0<a<1$$ is a minimal size of cells and $$u_0$$ represents an initial density function, which is defined on the initial set illustrated by the blue segment in Figure 5. Since Diekmann et al. (1984) interpret the functional term $$4b(2x)u(t,2x)$$ in (8) as zero whenever $$x\geq \frac{1}{2}\ ,$$ $$u(t_2,x_2)$$ has to be determined from (8) with $$4b(2x_2)u(t_2,2x_2)=0\ .$$ However, the value $$u(t_1,x_1)$$ has to be determined from (8) with $$4b(2x_1)u(t_1,2x_1)\ne 0$$ (because $$x_1 < \frac{1}{2}$$) and $$u(t_1,2x_1)$$ has to be determined like $$u(t_2,x_2)\ .$$

## Numerical Methods for Delay Partial Differential Equations

The most often applied numerical techniques for delay partial differential equations are composed of two consecutive steps:

• discretization in the variable $$x\ ,$$
• integration in $$t\ .$$

In the first step, the partial derivatives with respect to $$x$$ are replaced by some approximations. For example, application of the finite difference method to (3) means replacing the partial derivative with respect to $$x$$ by the approximating operator $\frac{u(t,x_{i-1})-2u(t,x_i)+u(t,x_{i+1})}{h^2} \approx \frac{\partial^2 }{\partial x^2} u(t,x_i).$ Here, $$h$$ is a step-size in $$x$$-direction and $$x_i$$ are grid-points defined by $x_i=ih, \quad i=0,1,\dots ,N, \quad h=\frac{\pi}{N}.$ This discretization in $$x$$ results in the following system of ordinary delay differential equations: $\tag{16} \frac{\partial u(t,x_i)}{\partial t} =D \frac{u(t,x_{i-1})-2u(t,x_i)+u(t,x_{i+1})}{h^2} + r u(t,x_i) \Bigg(1-\int_{-\tau }^0u(t+s,x_i)d\eta (s)\Bigg),$

where $$u(t,x_0)$$ and $$u(t,x_N)$$ are known from the boundary conditions in (11). Since, at this stage, $$t$$ is a continuous variable (not yet discretized), systems like (16) are called semi-discrete systems, and the first step is called the process of semi-discretization, or method of lines.

The finite difference methods (e.g. illustrated by (16)) are the most often applied numerical methods for the process of semi-discretization of delay partial differential equations, see e.g. van der Houwen et al (1986), Higham & Sardar (1995), and Zubik-Kowal & Vandewalle (1999). However, the Galerkin finite elements method has been successfully applied to (2) by Rey & Mackey (1993). Moreover, recent papers by Zubik-Kowal (2000), Mead & Zubik-Kowal (2005), and Jackiewicz & Zubik-Kowal (2006) show that pseudospectral methods solve delay partial differential equations with exponential accuracy; that is, their errors decay at exponential rates. Therefore, the accuracies of the finite difference and finite elements methods do not even come close to this exponential accuracy. Pseudospectral methods for non-delay and non-functional partial differential equations are investigated e.g. by Canuto et al. (1988).

In the second step of solving delay partial differential equations, the resulting semi-discrete systems are integrated in $$t\ .$$ These systems are composed of stiff ordinary delay differential equations. Numerical methods for solving such systems have been deeply investigated by Bellen & Zennaro (2003) (see also Stiff Delay Equations). The monograph by Bellen & Zennaro (2003) provides also deep investigations of numerical methods for solving ordinary delay differential equations, which are not necessarily results from semi-discretization of delay partial differential equations (see Delay Differential Equations).

Numerical methods for delay partial differential equations bring specific difficulties, which do not appear for equations without delays. For example, because of the delays in (1), previously computed approximations to $$u(\eta ,\xi)\ ,$$ for all $$\xi$$ and $$\bar{t}-\max\{ \tau_1,\tau_2\}\leq \eta \leq \bar{t}\ ,$$ have to be stored so that the next approximations for $$t\geq \bar{t}$$ can be computed. Since the transformed $$x$$-variable in (2) may not be included in the $$x$$-mesh, additional approximations have to be constructed for $$u(t,\alpha x)\ .$$ The integrals in (3) have to be approximated for all grid-points in the $$(t,x)$$-domain. The infinite domain of the integration in (5) has to be decomposed by using the structure of the kernel $$\gamma$$ so that the infinite integrals can be approximated by integrals over finite intervals. More details about numerical solutions for (4)-(5) are described by Zubik-Kowal (2006).