Diffusion on a Hypersphere: Application to the Wright-Fisher model
aa r X i v : . [ m a t h . P R ] D ec Diffusion on a Hypersphere: Application to theWright-Fisher model
Kishiko MaruyamaYoshiaki Itoh
The Institute of Statistical Mathematics and The Graduate University for AdvancedStudies, 10-3 Midori-cho, Tachikawa, Tokyo 190-8562, Japan
Abstract.
The eigenfunction expansion by Gegenbauer polynomials for the diffusionon a hypersphere is transformed into the diffusion for the Wright-Fisher model witha particular mutation rate. We use the Ito calculus considering stochastic differentialequations. The expansion gives a simple interpretation of the Griffiths eigenfunctionexpansion for the Wright-Fisher model. Our representation is useful to simulate theWright-Fisher model as well as Brownian motion on a hypersphere.PACS numbers: 02.50.Ey, 87.10.Ca iffusion on a Hypersphere: Application to the Wright-Fisher model
1. Introduction
Here we show that the diffusion on a hypersphere [3] is transformed into the diffusion forthe Wright-Fisher model with a particular mutation rate [7, 8, 9], by using the relation, x i = y i , where x i ’s denote the relative abundance of alleles and y i ’s denote the positionof a particle of the diffusion on a hypersphere. Diffusion on a sphere has been appliedto various problems in physics, chemistry, mathematics, etc. [4, 22, 28, 1]. The Wright-Fisher model is widely used to study the random sampling effect in population genetics[5, 6, 16, 26].We consider stochastic differential equations in the Ito sense [20]. We represent theWright-Fisher model for k alleles as well as the diffusion on ( k − (cid:16) k (cid:17) independent standard Brownian motion b ij ( t ) i > j with the skewsymmetry b ij ( t ) + b ji ( t ) = 0 [11, 13, 14, 18]. Our representation is useful to simulatethe Wright-Fisher model [19, 23]. There are effective simulation methods to simulate aBrownian motion on a hypersphere for example [2]. Our representation is simple andmay be useful to simulate the Brownian motion on a hypersphere.
2. Stochastic differential equation for the Wright-Fisher model
In the Wright-Fisher model each of the genes of the next generation is obtained by arandom choice among the genes of the previous generation and that the whole populationchanges all at once. In Moran’s model [21] it is supposed that there are 2 N individualseach formed from k alleles A , A , ..., A k , and that each instant at which the state of themodel may change, one individual of the alleles, chosen at random, dies and is replacedby a new individual which is A i with probability k i N , where k i is the abundance ofthe allele A i . It is supposed that the probability of any individual ”dying” during aninterval ( t, t + dt ) and then being replaced by a new individual is λdt . Hence the meannumber of such events in unit time is λ and the mean length of a generation is λ − . Thefollowing model of interacting particles [12] is equivalent with the Moran model.Consider a population of N particles each of which is one of k types, A , A , ..., A k .The types may represent species, alleles, genotypes or other classification. We thenconsider random interactions between particles, which are assumed to occur at the rate λdt per time interval ( t, t + dt ) for each particle. If a pair of particles of different types i and j interact, then after the interaction the both particles are the type i with probability1/2 and the type j with probability 1/2. If the type of the interacting particles are thesame, no change occurs.We can approximate the behavior of our interacting particle system by a stochasticdifferential equation (1). In it, the relative abundance of type i increases by c q x i ( t ) x j ( t ) db ij ( t ) and decreases by c q x i ( t ) x j ( t ) db ij ( t ) with the interaction of theparticles of type j , where c = q λ/ N . Hence our interacting particle systemautomatically makes the following equation (1), which has the genetic drift matrix withthe elements c x i ( t )( δ ij − x j ( t )) dt for i, j = 1 , , ..., k as covariances. iffusion on a Hypersphere: Application to the Wright-Fisher model i, j = 1 , , ..., k , consider dx i ( t ) = k X j =1 ,i = j c q x i ( t ) x j ( t ) db ij ( t ) (1)with b ij ( t ) + b ji ( t ) = 0, where b ij ( t ) ( i > j ) are mutually independent one-dimensionalBrownian motion with the mean 0 and the variance t [11, 13, 14, 18]. This representationof the Wright-Fisher model [11] is applied to make computer simulations on models inpopulation genetics [19, 23]. Here we make use of this equation (1) to discuss the exactsolutions of the diffusion equations considering the diffusion on a unit sphere.
3. Diffusion on hypersphere
The isotropic diffusion on the ( k − dy i ( t ) = − c k − y i ( t ) dt + c k X j =1 y j ( t ) db ij ( t ) (2)with b ij ( t ) + b ji ( t ) = 0, for i = 1 , , ..., k , where b ij ( t ) ( i > j ) are mutually independentone-dimensional Brownian motion with the mean 0 and the variance t . Let dy =( dy , ..., dy k ). By using the Ito calculus [20], we have for the stochastic differential d ( X i y i ( t ) ) = 0 . (3)Hence starting from a point on the unit hypersphere, the trajectory is on the unithypersphere P i y i ( t ) = 1. Consider the dot product for a unit vector l = ( l , l , ..., l k ), l · dy = X i l i [ c X j ( y j ( t ) db ij ( t )) + − c k − y i ( t ) dt ] . (4)On the tangent hyperplane at y = ( y , y , ..., y k ) of the hypersphere, for l which isperpendicular to y , we have l · dy = c X i l i ( X j y j ( t ) db ij ( t )) . (5)We have the expectation E ( l · dy ) = 0 . (6)We have V ar ( l · dy ) = ( c X i>j ( l i y j − l j y i ) dt = ( c dt, (7)since X i>j ( l i y j − y i l j ) = X i>j [( l i y j ) + ( y i l j ) − l i y j y i l j ] + ( X i l i y i ) = X i>j [( l i y j ) + ( y i l j ) ] + X i ( l i y i ) = ( X i l i )( X i y i ) = 1 . (8) iffusion on a Hypersphere: Application to the Wright-Fisher model V ar ( l · dy ) does not depend on the direction l at the tangential point y of the tangentialplane. Hence we see the diffusion is isotropic and we have the Fokker-Planck equationon the hypersphere for the Laplace-Beltrami operator [3] in space S k − , ∂∂t ρ ( y , t | y ′ ,
0) = D △ S k − ρ ( y , t | y ′ ,
0) (9)where D = c △ S k − is for k − θ , θ , . . . , θ k − of the unit vector y = ( y , ..., y k ) of S k − defined by the relations y = sin θ k − ... sin θ sin θ (10) y = sin θ k − ... sin θ cos θ (11)... (12) y k − = sin θ k − cos θ k − (13) y k = cos θ k − , (14)where 0 ≤ θ < π and 0 ≤ θ i < π for i = 1. The integration measure in S k − is definedas dy = A k − sin k − θ k − ... sin θ dθ ...dθ k − , where A k − = 2 π k/ / Γ( k/
2) is the surfaceof the sphere S k − . With this normalization R S k − dy = 1.The solution [3], the transition probability density for the isotropic diffusion onhypersphere, is obtained by using the definitions and notations of the chapter IX ofthe book [24]. The Gegenbauer polynomial C pL is defined as the coefficient of h L in thepower-series expansion of the function(1 − th + h ) − p = ∞ X L =0 C pL ( t ) h L . (15)The transition probability density from y ′ (with n − θ ′ , . . . , θ ′ n − )at 0 to y (with θ , . . . , θ n − ) at t > ρ ( y , t | y ′ , A k − ∞ X L =0 L + k − k − C k/ − L ( y · y ′ ) exp( − DL ( L + k − t ) , (16)for the dot product y · y ′ = P i y i y ′ i , where C ( k/ − L ( z ) = ⌊ L/ ⌋ X j =0 ( − j Γ( L − j + k/ − k/ − j !( L − j )! (2 z ) L − j . (17)
4. Solution for the Wright-Fisher model with parent-independent mutation
Let A , ..., A k denote k allele types in a population. The general neutral alleles modelhas a probability u ij of a mutation from A i to A j . An exact solution is known for thecase of the 0 mutation rate case in three dimensions [16]. Assuming parent-independentmutation [7], u ij = u j . an expansion of the transition probability density [7], from the iffusion on a Hypersphere: Application to the Wright-Fisher model x ′ = ( x ′ , ..., x ′ k ) to the relative abundances x = ( x , ..., x k ), isobtained for the following Fokker-Planck equation (18), ∂p ( x , t | x ′ , ∂t = k − X i =1 ∂∂x i ( 12 k − X i =1 ∂x i ( δ ij − x j ) ∂x j − M i ) p ( x , t | x ′ ,
0) (18)where M i = ε i − µx i , µ = P ki =1 ε i . For ε i > , i = 1 , , ..., k , the stationary density[27, 25] is given byΓ( µ ) k Y i =1 x ε i − i / Γ( ε i ) . (19)The solution of equation (18) is given by an expansion in orthogonal polynomials [7].For the case ε i = ε , i = 1 , ..., k , we have the solution [7], p ( x , t | x ′ , µ ) k Y i =1 x ε − i Γ( ε ) ( ∞ X n =0 exp {− n ( n − t − µnt } Q n ( x , x ′ )) (20)for x j > , j = 1 , , ..., k , where Q n is given by Q n ( x , x ′ ) = ( µ + 2 n − n !) − n X m =0 ( − n − m nm ! ( µ + m ) ( n − ξ m (21)with ξ m = µ ( m ) Γ( ε ) k X l + ... + l k = m, ≤ l j ,j =1 ,...,k m ! l ! ..., l k ! Π kj =1 ( x j x ′ j ) l j Γ( l j + ε ) , (22)for the notation α ( m ) = α ( α + 1) · · · ( α + m − Q ( x , x ′ ) = 1 , (23) Q ( x , x ′ ) = ( µ + 1)( ξ − , (24) Q ( x , x ′ ) = 12 ( µ + 3)[( µ + 2) ξ − µ + 1) ξ + µ )] , (25)with ξ = µ k X j =1 x j x ′ j ε (26) ξ = µ ( µ + 1) k X j =1 ( ( x j x ′ j ) ε ( ε + 1) + 2 X i 0) for equation (18) with ε i = 1 / i = , ..., k on the relative abundances x = ( x , ..., x k ) starting from x ′ = ( x ′ ...., x ′ k ) isrepresented by using the transition probability density ρ ( y , t | y ′ , y = ( y , ..., y k ) with spherical coordinates θ , ..., θ k of ( k − y ′ = ( y ′ , ..., y ′ k ) with spherical coordinates θ ′ , ..., θ ′ k , as p ( x , t | x ′ , (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∂ ( x , ..., x k − ) ∂ ( θ , ..., θ k − ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) dθ ...dθ k − (30)= X y ′ i = x ′ i , ≤ y ′ i ,y i = x i ,i =1 ,...,k ρ ( y , t | y ′ , (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∂ ( y , ..., y k − ) ∂ ( θ , ..., θ k − ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) dθ ...dθ k − . (31)We have (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∂ ( x , ..., x n − ) ∂ ( θ , ..., θ n − ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 2 k − ( k Y i =1 y i ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∂ ( y , ..., y k − ) ∂ ( θ , ..., θ k − ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (32)Considering equation (16), ( ?? ), (17), (32) we have p ( x , t | x ′ , k/ π k/ ( k Y i =1 x − / i ) 2 − k X y ′ i = x ′ i , ≤ y ′ i ,y i = x i ,i =1 ,...,k ∞ X L =0 L + k − k − C k/ − L ( y · y ′ ) exp( − DL ( L + k − t ) . (33)Thus we proved the solution (33) is equivalent with the solution (20) for the case ε = 1 / 2. Let us observe the two forms of the solution. Consider the eigen values of thetwo diffusion equations (9) and (18). The Caillol solution [3] of (9), expanded by theGegenbauer polynomials, are applied in equation (33). The Gegenbauer polynomials forodd L in equation (33) cancel with each other as we can see from equation (17). Hencethere is no term for odd L in equation (33), [17]. Put L = 2 n , and take D = 1 / c = 1 in(7)), for particular mutation rate ε i = 1 / , i = 1 , ..., k . Then the exponential functionsexp( − DL ( L + k − t ) in equation (33) coincide with those of Griffiths solution (20).The above Q ( x , x ′ ), Q ( x , x ′ ), and Q ( x , x ′ ) in the solution (20) are obtained from theterms in the expansion (33) by using x i = y i and x ′ i = y ′ i for i = 1 , ..., k , as we cancalculate easily. Also we seeΓ( k/ π k/ ( k Y i =1 x − / i ) = Γ( µ ) k Y i =1 x ε − i Γ( ε ) . (34)Our study naturally connects the diffusion on a hypersphere with the Wright-Fishermodel for a particular mutation rate of the parent-independent mutation.Remark 1. The connection between diffusion on a hypersphere and the 1-dimensional Wright-Fisher diffusion has been derived before (p338 in [15] ). Theequation (13.33) in [15] corresponds to equation (16). They derive the spectral expansion iffusion on a Hypersphere: Application to the Wright-Fisher model Acknowledgement. The authors thank Hidenori Tachida for the references inpopulation genetics. The authors thank Shinichi Kotani, Tohru Ogawa and Geoffrey A.Watterson for helpful discussion and suggestion on the symmetry of special functions.The authors thank the reviewers for helpful comments and suggestions. [1] Brillinger, D. R. (1997). A Particle Migrating Randomly on a Sphere. Journal of TheoreticalProbability, 10(2), 429-443.[2] Carlsson, T., Ekholm, T., & Elvingson, C. (2010). Algorithm for generating a Brownian motion ona sphere. Journal of physics A: Mathematical and theoretical, 43(50), 505001.[3] Caillol, J. M. (2004). Random walks on hyperspheres of arbitrary dimensions. Journal of PhysicsA: Mathematical and General, 37(9), 3077.[4] Debye P. (1929) Polar molecules, Dover, London.[5] Ewens, W. J. (2004) Mathematical Population Genetics 1: I. Theoretical Introduction.SpringerScience & Business Media, Vol. 27. .[6] Fisher, R. A. (1922) On the dominance ratio, Proceedings of the Royal Society of Edinburgh, 42,321-341.[7] Griffiths, R. C. (1979). A transition density expansion for a multi-allele diffusion model. Advancesin Applied Probability, 310-325.[8] Griffiths, R. C. (1980). Allele frequencies in multidimensional Wright-Fisher models with a generalsymmetric mutation structure. Theoretical Population Biology, 17(1), 51-70.[9] Griffiths, B. (2009). Stochastic processes with orthogonal polynomial eigenfunctions. Journal ofcomputational and applied mathematics, 233(3), 739-744.[10] Griffiths, R. C., and Spano, D. (2013). Orthogonal polynomial kernels and canonical correlationsfor Dirichlet measures. Bernoulli, 19(2), 548-598.[11] Itoh, Y. (1979). Random collision process on oriented graph. Institute of Statistical Mathematics(Japan), Research Memorandum, (154), 1-20.[12] Itoh, Y. (1979). Random collision models in oriented graphs. Journal of Applied Probability, 36-44.[13] Itoh, Y. (1993). Stochastic model of an integrable nonlinear system. Journal of the Physical Societyof Japan, 62(6), 1826-1828.[14] Itoh, Y., Mallows, C., and Shepp, L. (1998) Explicit sufficient invariants for an interacting particlesystem.” Journal of Applied Probability, 35, 633-641.[15] Karlin, S. and Taylor, H. (1981) A second course in stochastic processes, Academic Press, NewYork.[16] Kimura, M. (1956). Random genetic drift in a tri-allelic locus; exact solution with a continuousmodel. Biometrics, 12(1), 57-66.[17] Maruyama, K. (1992) Generalized Random Frequency Modulation Model and Its Applications,PhD Dissertation, The Graduate University for Advanced Studies (in Japanese).[18] Maruyama, K. and Itoh, Y. (1991) Stochastic differential equations for the Wright-Fisher model,Proceedings of the Institute of Statistical Mathematics, 39, 47-52 (in Japnese with English sumary).[19] Maruyama, T., and Nei, M. (1981). Genetic variability maintained by mutation and overdominant iffusion on a Hypersphere: Application to the Wright-Fisher model8