Prepared using rsifpublic.cls

J. R. Soc. Lond. Interface (0000) DOI: 10.1098/rsif0000.0000

Supplementary material to Mathematical models of the VEGF receptor and its role in cancer therapy Tom´ as Alarc´ on† and Karen M. Page Bioinformatics Unit, Department of Computer Science, University College London, London WC1E 6BT, UK These notes, accompanying the paper Mathematical models of the VEGF receptor and its role in cancer therapy, contain the proof of Equations (10) and (12), which constitutes a generalisation of the proof presented originally by Kubo et al, J. Stat. Phys. 9, 51-96 (1973) for one-dimensional systems. This proof is also included with a pedagogic purpose, as this method, well-known among the statistical physics community, is relatively less known among applied mathematicians. For completeness, we also include an alternative method for the derivation of Eqs. (10) and (12) based on notions of stochastic calculus and stochastic differential equations.

Hereafter, when we refer to Eqs. (10) and (12), we mean Eqs. (10) and (12) of the main body of the paper Mathematical models of the VEGF receptor and its role in cancer therapy, unless stated otherwise.

of the solution of the ME. In this paper we will restrict ourselves to the system of equations for the first moments (q1 (t)) and for the elements of the covariance matrix (q2 (t)), which will enable us to study the average behaviour of the system and the fluctuations around it. The details of the derivation of the equation for Q(~u, t) are rather technical and it is given in full detail in Section 1.1. The evolution equation for Q(~u, t) is given by

1. WKB method The starting point for our derivation of Eqs. (10) and (12) of the paper that these notes accompany is the Master Equation (ME) written in WKB form: X ~ r ∂ 1 ∂ψ(~x, t) (e− N · ∂~x − 1)a(~x, ~r, t)ψ(~x, t), = N ∂t

1 ∂Q(~u, t) = N ∂t Z ∞ 1 X  −i u~ ·~r N e −1 d~v Q(~u − ~v , t)w(~v , ~r, t), (2π)d −∞

(1)

~ r

~ r

∂

where we have used that e−~r· ∂~x is the generator of the translations in the space of states of the system. Let us consider the cumulant-generating function, given by q(~u, t) = ln Q(~u, t), where Q(~u, t) is the characteristic function of ψ(~x, t), defined as its d-dimensional Fourier transform: Z ∞ Q(~u, t) = d~x ei~u·~x ψ(~x, t). (2)

(4) where w(~v , ~r, t) is the Fourier transform of a(~x, ~r, t). Eq. (4) will serve as the basis for obtaining the equations for the cumulants, which is our main aim. 1.1. Derivation of Eq. (4) To obtain the equation for Q(~u, t), we take the Fourier transform of the ME (1):

−∞

The (multivariate) cumulants, qn (t), can be obtained from q(~u, t) as the coefficients of the expansion: q(~u, t) =

∞ n X i n ~u qn (t), n! n=1

XZ ∞ ~ r ∂ 1 ∂Q(~u, t) = d~x ei~u·~x (e− N · ∂~x − 1)a(~x, ~r, t)ψ(~x, t). N ∂t −∞ ~ r (5) The integral on the right hand side of Eq. (5) is:

(3)

where ~un stands for the n-adic product defined as (~un )i1 i2 ..in ≡ ui1 ui1 ..uin and the symbol“ ” denotes full contraction. The procedure we follow next, which is based on the work of Kubo et al (1973), is to write down an equation for Q(~u, t) and construct asymptotic approximations for Q(~u, t), q(~u, t) and, finally, qn (t). These approximations will provide systems of ODEs for the cumulants, and, consequently, for the moments

Z

∞

~ r

∂

d~x ei~u·~x (e− N · ∂~x − 1)a(~x, ~r, t)ψ(~x, t) =

−∞ ∞ X n=0

Z

∞

d~x e −∞

i~ u·~ x



~r ∂ − · N ∂~x

n a(~x, ~r, t)ψ(~x, t). (6)

1

Article submitted to The Royal Society

2

Modelling the VEGF receptor: Supplementary material. T. Alarc´on and K.M. Page where q(~u − ~v , t) has been replaced by q(~u − ~v , t) = e~u·∂~v q(−~v , t) and q (n) (−~v , t) ≡ (∂~u )n q(~u, t)|~u=−~v . Using Eq. (11), Eq. (4) reads:

In order to proceed further we use the identity:

ei~u·~x 

 −

~r ∂ · N ∂~x

n a(~x, ~r, t)ψ(~x, t) =

~r · ~u ~r ∂ −i − · N ∂~x N

n

ei~u·~x a(~x, ~r, t)ψ(~x, t), (7)

which is proved in Section 1.3. Using Eq. (7) and the binomial formula Eq. (6) reads: ∞

Z

~ r

∂

d~x ei~u·~x (e− N · ∂~x − 1)a(~x, ~r, t)ψ(~x, t)

−∞ ∞ Z X

∞

n X

n! k!(n − k)! n=1 −∞ k=0 k  n−k !  ~r · ~u ~r ∂ ei~u·~x a(~x, ~r, t)ψ(~x, t). −i − · N ∂~x N

=

d~x

(8) For any value of n, the corresponding sum over k yields only one term without derivatives (corresponding to k = 0). All the other terms (0 < k ≤ n) consist of terms that, upon integration by parts, produce boundary terms that vanish under appropriate conditions. Thus, for any n, there is only one term in the corresponding k-sum. Taking this into account, we obtain:

∞ 1 X in n ~u q˙n (t) = N n=1 n! Z ∞ 1 X  −i u~ ·~r Q(~u − ~v , t) N e −1 d~v w(~v , ~r, t) (2π)d Q(~u, t) −∞ ~ r Z  h u~ ·~r 1 X ∞ −i N = − 1 d~ v e (2π)d −∞ ~ r    ∞ j X i ~uj hj (−~v , t) eq(−~v,t) w(~v , ~r, t) , exp  j! j=1

(12) where, for convenience, we have defined the fol-
lowing quantity hj (−~v , t) ≡ i−j q (j) (−~v , t) − q (j) (0, t) and we have used q (0) (−~v , t) = q(−~v , t). Let us focus now on the first cumulant q1 (t) = hx(t)i, which is given by the O(~u) terms in Eq. (12) (of the supplementary materials). Expanding the exponentials within the right hand side of Eq. (12) (of the supplementary materials) and keeping only first order terms we obtain: Z ∞  X ~ r d~v eq(−~v,t) w(~v , ~r, t)m0 (, ~v , t), (2π)d −∞ ~ r (13) where  ≡ N −1 and mk (, ~v , t) is defined by: q˙1 (t) =

Z

∞

~ r

∂

d~x ei~u·~x (e− N · ∂~x − 1)a(~x, ~r, t)ψ(~x, t) =

−∞

 n Z ∞ ∞ X 1 ~r · ~u −i d~xei~u·~x a(~x, ~r, t)ψ(~x, t). n! N −∞ n=1 (9)



 ∞ j ∞ n X X i j i n exp  ~u hj (−~v , t) = ~u mn (, ~v , t). j! n! j=1 k=0

Substituting Eq. (9) in Eq. (5), rearranging terms, and recalling that the Fourier transform of the product of two functions equals the convolution of the corresponding Fourier transforms, we finally obtain: 1 ∂Q(~u, t) = N ∂t Z ∞ 1 X  −i u~ ·~r N d~v Q(~u − ~v , t)w(~v , ~r, t) e −1 (2π)d −∞ ~ r

(10) where w(~v , ~r, t) is the Fourier transform of a(~x, ~r, t). 1.2. Derivation of Eqs. (10) and (12) Before proceeding forward, we note that:

Q(~u − ~v , t) = e

q(~ u−~ v ,t)

= exp

Prepared using rsifpublic.cls

! ∞ n X i n (n) ~u q (−~v , t) , n! n=0 (11)

(14) Eq. (13) needs to be properly balanced with respect to the small parameter . Such balance is achieved if, to leading order, m0 (, ~v , t) = O(0 ), which is, in fact, satisfied as m0 (, ~v , t) = 1 (see Eq. (14)), and q1 (t) = O(0 ). Likewise, for n = 2 we obtain:  q˙2 (t) = 2

Z 1 X ∞ d~v ( m1 (, ~v , t)~r+ (2π)d −∞ ~ r  2 m0 (, ~v , t)~r2 eq(−~v,t) w(~v , ~r, t).(15) 2

Balancing of Eq. (15) leads to the leading-order scaling q2 (t) = O() and m1 (, ~v , t) = O(). These scaling relationships can be checked to be consistent with Eq. (14). In general, we can show that the scaling substitution qn (t) = n−1 qn1 (t) + n qn2 (t) + O(n+1 ) and mn (, ~v , t) = n mn1 (~v , t) + n+1 mn2 (~v , t) + O(n+2 ) is consistent with Eq. (14) and leads to balanced equations for the cumulants qn (t).

Modelling the VEGF receptor: Supplementary material. T. Alarc´on and K.M. Page To leading order, the equation for q1 (t) reads:

3

1.3. Proof of identity Eq. (7) In this Section we prove the identity Eq. (7):

Z ∞ 1 X q˙1 (t) = ~r d~v eq(−~v,t) w(~v , ~r, t), (2π)d −∞

(16)

~ r

ei~u·~x

whereas the corresponding equation for q2 (t) is given by:  q˙2 (t) = 2



Z 1 X ∞ d~v ( m1 (, ~v , t)~r (2π)d −∞ ~ r  2 2 q(−~v,t) w(~v , ~r, t). (17) + ~r e 2

∞ X

e

  in+j j (∂~u ) ~u(n+j) q(n+j) (t) |~u=−~v . (n + j)! n=1 (18) For j = 1, from Eq. (14), we have that m1 = h1 . Taking into account the scaling substitution for mn and qn , we obtain the following O() approximation:
i  m11 (~v , t) =  ∂~v ~v 2 q21 (t) . 2

(19)

Thus, substituting Eq. (19) into Eq. (17) and keeping only the leading order contributions, the equation for q21 (t) reads:

−

~r ∂ · N ∂~x

n

~r · ~u ~r ∂ −i − · N ∂~x N

a(~x, ~r, t)ψ(~x, t) = n

ei~u·~x a(~x, ~r, t)ψ(~x, t),

The proof will proceed by induction. Proof. Let n ∈ N. n = 1. Trivial. n > 1. Let us assume Eq. (7) holds for n − 1. Then:

The leading order approximation of Eq. (17) is obtained as follows. From the definitions of the quantities hj and mj we obtain the following expression for hj :

hj (−~v , t) = i−j



i~ u·~ x



~r ∂ − · N ∂~x

n a(~x, ~r, t)ψ(~x, t) =

!  n−1 ~r ∂ ~r ∂ i~ u·~ x a(~x, ~r, t)ψ(~x, t) + − · e − · N ∂~x N ∂~x    n−1 ~r ∂ ~r · ~u i~ u·~ x a(~x, ~r, t)ψ(~x, t) = e − · −i N N ∂~x !  n−1 ~r ∂ ~r · ~u ~r ∂ i~ u·~ x − · −i − · e a(~x, ~r, t)ψ(~x, t) + N ∂~x N ∂~x N   n−1 ~r · ~u ~r · ~u ~r ∂ −i −i − · ei~u·~x a(~x, ~r, t)ψ(~x, t) = N N ∂~x N  n−1  ~r ∂ ~r · ~u ~r · ~u ~r ∂ − · −i −i − · N ∂~x N N ∂~x N ei~u·~x a(~x, ~r, t)ψ(~x, t) = n  ~r · ~u ~r ∂ −i ei~u·~x a(~x, ~r, t)ψ(~x, t). QED. − · N ∂~x N

Z
1 X ∞ d~v i~v · q21 (t)~r + ~r2 eq(−~v,t) w(~v , ~r, t). d (2π) −∞ ~ r (20) 2. Alternative method Furthermore, the scaling substitution for qn (t) and Here we explain a second method for the derivation Eq. (3) lead to: of Eqs. (10) and (12) based on concepts of stochastic calculus and stochastic differential equations (SDEs). 1 Those familiar with these disciplines may find this (21) q(~u, t) = ρ(, ~u, t) second alternative more straightforward.  ∞ X ∞ n To proceed further, we follow Allen (2003), in particX
i n ρ(, ~u, t) = ~u n+k−1 qnk (t) .(22) ular her derivation of a set of Ito SDEs for interacting n! n=1 k=1 populations. As detailed in Sections 2 and 3 of the ~ defined as the paper, the state vector of the system, X, This almost completes our derivation of the evolution vector whose components are the number of particles equations for q1 (t) and q2 (t) to leading order. Eq. (22) of each of the (chemical) species, during the interval can be rewritten as ~ →X ~ + ~r, (t, t + ∆t) changes by an amount ~ri , i.e. X ~ with probability W (X, ~r, t)∆t. Thus the mean of the ~ = ~r, is given by: change in the state vector, ∆X 2 ρ(, ~u, t) = ρ1 (~u, t) + O( ) (23) X ~ = ∆t ~ ~r, t) ≡ µ ρ1 (~u, t) = i~u · q11 (t). (24) h∆Xi ~rW (X, ~ ∆t. (25) q˙21 (t) =

~ r

Substituting Eqs. (23) and (24) into Eq. (16) we obtain that, to leading order, the equations for q11 (t) and q21 (t) are Eqs. (10) and (12). Prepared using rsifpublic.cls

~ defined as Likewise, the covariance matrix of ∆X, T T 2 ~ = h∆X∆ ~ X ~ i−µ Cov(∆X) ~µ ~ ∆t , is given by:

4

Modelling the VEGF receptor: Supplementary material. T. Alarc´on and K.M. Page which, to leading in order in dt, gives: ~ = ∆t Cov(∆X)

X

~ ~r, t) ≡ C∆t, ~r~rT W (X,

(26)

~ r

where we have assumed ∆t → 0 and neglected higher order terms. The matrix C is symmetrical so has real eigenvalues. We will assume that, in addition, C is positive definite. This assumption allows us to define a matrix, B, such that C = BB T . ~ is suffiFurthermore, according to Allen (2003), if X ~ has a norciently large, it is possible to assume that ∆X mal distribution with mean µ ~ ∆t and covariance matrix ~ = N (~ BB T ∆t, i.e. ∆X µ∆t, BB T ∆t). Taking this into account we define a random vector ~η = N (0, I), √ where I is the identity matrix, such that µ ~ ∆t + B ∆t~η = ~ + ∆t) = X(t) ~ ~ can be N (~ µ∆t, BB T ∆t). Thus X(t + ∆X written as: √ ~ + ∆t) = X(t) ~ X(t +µ ~ ∆t + ∆tB · ~η . (27) √ Taking now the limit dt ≡ ∆t → 0 and assuming that ∆tηi → dWi , where Wi is the Wiener process, the system converges (in the mean square sense) to the system of Ito SDEs: ~ ~ dX(t) =µ ~ dt + B · dW.

(28)

Recalling the homogeneity property of the transition ~ ~r, t) = N a(~x, ~r, t), with ~x = X/N ~ rates W (X, , Eq. (28) reads: d~x(t) = dt

X

~ra(~x, ~r, t) +

~ r

1 ~ BdW. N

(29)

Since B/N = O(N −1/2 ) in the limit N → ∞ the noise term in Eq. (29) can be neglected. In this case the process is almost deterministic, i.e. ~x ' h~xi. Under these two assumptions, Eq. (10) follows from Eq. (29): dh~xi X = ~ra(h~xi, ~r, t). dt

(30)

~ r

Eq. (12) can be obtained from Eq. (28) by applying the multidimensional Ito formula (Oksendal, 2003). The (multidimensional) Ito formula, which can be viewed as the stochastic generalisation of the chain rule in ~ is described by the standard calculus, states that if X ~ is also an Ito Ito equation Eq. (28), then Y = f (t, X) process given by

dY =

2 ∂f ∂f ~ + 1 ∂ f : dXd ~ X ~ T, · dX dt + ~ ~ X ~ ∂t 2 ∂ X∂ ∂X

(31)

where the colon denotes double contraction. In the ~X ~ T − hXih ~ Xi ~ T or, in components, present case Y = X Yij = Xi Xj − hXi ihXj i. After some algebra it is easy to check that Eq. (31) reduces to

dYij

1 = Xi dXj + Xj dXi + (dXi dXj + dXj dXi ) 2 −hXi iµj dt − hXj iµi dt, (32)

Prepared using rsifpublic.cls

dYij = (Xi − hXi i)µj dt + (Xj − hXj i)µi dt X + (µi Bjk + µj Bik ) dtdWk l

1X + (Bik Bjl + Bjk Bil )dWk dWl 2 k,l X + (Xi Bjk + Xj Bik )dWk .

(33)

k

To proceed further, we take the Taylor expansion ~ keeping the of the functions µi , i = 1, 2, ... around hXi, terms up to first order: ~ ~r, t) ' µi (hXi, ~ ~r, t) + µi (X,

X ∂µi ~ ~ (Xk − hXk i) ∂Xk X=hXi k

Introducing this Taylor expansion in Eq. (33) and taking the average, we obtain: dhYij i = dt

X k

+

X

 ∂µi ∂µj hYik i ~ ~ + hYjk i ~ ~ ∂Xk X=hXi ∂Xk X=hXi

h(µi Bjk + µj Bik ) dtdWk i

l

1X h(Bik Bjl + Bjk Bil )dWk dWl i 2 k,l X + h(Xi Bjk + Xj Bik )dWk i. +

(34)

k

To carry on with our derivation we assume statisti~ and dW, ~ i.e. hXi dWj i = cal independence between X hXi ihdWj i. That this is, in fact, the case can be argued as follows. Consider Eq. (27), which constitutes the discretisation of the Ito equation Eq. (28)√according to the Euler method. Moreover hXi dWj i ' ∆thXi ηj i for ∆t → 0. Since ηj = N (0, 1), ηj is independent of ~ Xi . Furthermore, according to Eq. (27), whilst X(t) depends on the randomly generated series of values of ~η (i.e. on the realisation of the Wiener process being generated) up to time t − ∆t, is independent of the value of ~η generated at time t. Thus, we can write hX j i which in the limit ∆t → 0 √i ηj i = hXi ihη√ implies ∆thXi ηj i ' ∆thηj ihXi i = hXi ihdWj i. This ~ result can be easily extended to a function of X(t) and ~ ~ ~ t, f (X(t), t): hf (X(t), t)dWj i = hf (X(t), t)ihdWj i. By an argument similar to the one explained in the previous paragraph we can generalise this result ~ ~ to write: hf (X(t), t)dWi dWj i = hf (X(t), t)ihdWi dWj i. These results, along with the properties of the Wiener process: hdWj i = 0 and hdWi dWj i = δij dt, allow us to write Eq. (34) as:   ∂µi dhYij i X ∂µj = hYik i ~ ~ + hYjk i ~ ~ dt ∂Xk X=hXi ∂Xk X=hXi k X ~ ~r, t)i. + ri rj hW (X, (35) ~ r

Modelling the VEGF receptor: Supplementary material. T. Alarc´on and K.M. Page Likewise, carrying out the corresponding Taylor ~ ~r, t) around hXi, ~ we expansion of the function W (X, obtain:   dhYij i X ∂µj ∂µi = hYik i + hY i X=h jk ~ ~ ~ ~ Xi dt ∂Xk ∂Xk X=hXi k   X X ∂W ~ ~r, t) + + ri rj W (hXi, hYkl i ~ ~  . ∂Xk ∂Xl X=hXi ~ r

k,l

(36) Using the scaling laws Yij = O(N ), µi = O(N ) and Xi = O(N ), we determine which terms of Eq. (36) vanish in the limit N >> 1. All the terms Eq. (36) are O(N ), except the term in the second derivatives of W which is O(1). Thus, to leading order, Eq. (36) reads:   ∂µj ∂µi dhYij i X = hYik i + hY i X=h jk ~ ~ ~ ~ Xi dt ∂Xk ∂Xk X=hXi k X ~ ~r, t). + ri rj W (hXi, (37) ~ r

which is Eq. (12). REFERENCES L.J.S. Allen. (2003). An introduction to stochastic processes with applications to biology Pearson Education, Upper Saddle River, NJ, USA. R. Kubo, K. Matsuo, and K. Kitahara. (1973). Fluctuation and relaxation of macrovariables J. Stat. Phys. 9, 51-96. B. Oksendal. (2003). Stochastic differential equations. An introduction with applications Springer-Verlag, Berlin, Germany.

Prepared using rsifpublic.cls

5