# Optimal Stochastic Control Problem for General Linear Dynamical Systems in Neuroscience

Authors: Yan Chen; Yingchun Deng; Shengjie Yue; Chao Deng
Journal: Advances in Mathematical Physics
 (2017)

---

## Abstract

This paper considers ad-dimensional stochastic optimization problem in neuroscience. Suppose the arm’s movement trajectory is modeled by high-order linear stochastic differential dynamic system in d-dimensional space, the optimal trajectory, velocity, and variance are explicitly obtained by using stochastic control method, which allows us to analytically establish exact relationships between various quantities. Moreover, the optimal trajectory is almost a straight line for a reaching movement; the optimal velocity bell-shaped and the optimal variance are consistent with the experimental Fitts law; that is, the longer the time of a reaching movement, the higher the accuracy of arriving at the target position, and the results can be directly applied to designing a reaching movement performed by a robotic arm in a more general environment.


---

## Body

## 1. Introduction

The effective control of neuronal activity is one of the most exciting topics in theoretical neuroscience, with great potential for applications in healthcare. Nowadays, the application of stochastic control methods in neuroscience is becoming a significant portion of the mainstream research. Among many researches, for example, we refer to Holden (1976) for the models of the stochastic activity of neural aggregates, Iolov et al. [1] with respect to the optimal control of single neuron spike trains, and Roberts et al. [2] with respect to the review of the application of the stochastic models of brain activity.In this paper, we study trajectory planning and control in human arm movements. When a hand is moved to a target, the central nervous system must select one specific trajectory among an infinite number of possible trajectories that lead to the target position. The content of this paper includes two parts: the first part is modeling the activities incorporating stochastic process, and the second part is quantifying task goals as cost functions and applying the sophisticated tools of optimal control theory to obtain the optimal behavior. Feng et al. [3] reviewed two optimal control problems at a different levels, neuronal activity control and movement control. They also derived the optimal signals for these two control problems. Li et al. [4] considered the robust control of human arm movements. Based on the fuzzy interpolation of an nonlinear stochastic arm system, they simplified the complex noise-tolerant robust control of the human arm tracking problem by solving a set of linear matrix inequalities using Newton’s iterative method via an interior point scheme for convex optimization. Singh et al. [5] modeled reaching movements in the presence of obstacles and solved a stochastic optimal control problem that consists of probabilistic collision avoidance constraints and a cost function that trades off between effort and end-state variance in the presence of a signal-dependent noise. For more details, we refer the reader to Campos and Calado [6], Berret et al. [7], and Mainprice et al. [8].Yet, all the above studies discussed 1-dimensional or lower-dimension space, and the neuronal activity or movement trajectory would be involved in a higher dimension space. In this paper, motivated by Feng et al. [3], we consider a stochastic control problem for arm movement within the framework of d-dimensional control space. Applying stochastic control theory, we solved the optimization problem explicitly and obtained the exact solution of the optimal trajectory, velocity, and the optimal variance.The remainder of this paper is organized as follows. Section2 introduces the basic model setup of the high-order linear stochastic dynamical systems for movement trajectory. In Section 3, we derive the explicit expressions for the optimal trajectory, velocity, and variance. In Section 4, we provided a 3-dimensional optimization example, and concluding remarks are given in Section 5.


## 2. Model Setup

### 2.1. The Integrate-and-Fire Model

In this subsection, we give the classical I & F (integrated-and-fire) model followed by Feng et al. [3], which describe the neuron activity.(1)dKt=-Kt-Krestνdt-dIsynt,t>0,with K(0)=Krest<Kthres and where ν is the decay time constant. The synaptic input current is (2)Isynt=a∑i=1pEit-b∑i=1qIit with Ei=Ei(t),t≥0 and Ij=Ij(t),t≥0 as Poisson processes with rates λE,i(t) and λI,j(t), a>0 and b>0 being the magnitude of each excitatory postsynaptic potential (EPSP) and inhibitory postsynaptic potential (IPSP); a cell receives excitatory postsynaptic potentials (EPSPs) at p synapses and inhibitory postsynaptic potentials (IPSPs) at q inhibitory synapses. Once K(t) crosses Kthre from below, a spike is generated and K is reset to Krest. This model is termed as the IF model.LetKrest=0,a=b,p=q, and use the usual approximation to approximate the IF models (see Feng et al. [3] and Zhang and Feng [9]); then (1) can rewriten as(3)dKt=-1νKtdt+aλtdt+aλαtdWt,where {W(t)}t≥0 is a standard Brownian motion, α>0 is a constant, and if α=0.5, it implies that the input is derived from a Poisson process. If α>0.5, it is the so-called the supra-Poisson inputs, and the other is the so-called sub-Poisson inputs if α<0.5. In addition, a larger α leads to more randomness for the synaptic inputs.


### 2.2. General Linear Stochastic Differential Equation

In this subsection, we extend the one-dimensional I & F model (3) to n-dimensional stochastic differential equations in which the solution process enters linearly. Such processes arise in estimation and control of linear systems, in systems, in economics, and in various other fields (see Liu [10]), as(4)dXt=AtXt+Λtdt+ΣtdWt,X0=ξ,where W is an r-dimensional Brown motion independent of the n-dimensional initial vector ξ, and n×n, n×1, and n×r matrices A(t), Λ(t), and Σ(t) are nonrandom, measurable, and locally bounded, respectively.Now we define ann×n matrix function Φ(t), satisfying the following matrix differential equation:(5)Φ˙t=AtΦt,Φ0=I,where I is the d×d identity matrix. We know that (3) has unique (absolutely continuous) solution defined for 0≤t<∞, and, for each t≥0, the matrix Φ(t) is nonsingular.By Itô’s rule, it is easily verified that(6)Xt=Φtξ+∫0tΦ-1sΛsds+∫0tΦ-1sΣsdWs;0≤t<∞solves (4).We suppose thatEξ2<∞ and introduce the mean vector and covariance matrix functions as follows:(7)mt=EXt,ρs,t=EXt-msXt-mtT,Vt=ρt,t.From (4), we can show that(8)mt=Φtm0+∫0tΦ-1sΛsds,(9)ρs,t=ΦsV0+∫0s∧tΦ-1uΣuΣuΦ-1u′duΦ′t hold for every 0≤s,t<∞. In particular, m(t) and V(t) solve the linear equations:(10)m˙t=Atmt+Λt,V˙t=AtΛt+ΛtA′t+ΣtΣ′t.


## 2.1. The Integrate-and-Fire Model

In this subsection, we give the classical I & F (integrated-and-fire) model followed by Feng et al. [3], which describe the neuron activity.(1)dKt=-Kt-Krestνdt-dIsynt,t>0,with K(0)=Krest<Kthres and where ν is the decay time constant. The synaptic input current is (2)Isynt=a∑i=1pEit-b∑i=1qIit with Ei=Ei(t),t≥0 and Ij=Ij(t),t≥0 as Poisson processes with rates λE,i(t) and λI,j(t), a>0 and b>0 being the magnitude of each excitatory postsynaptic potential (EPSP) and inhibitory postsynaptic potential (IPSP); a cell receives excitatory postsynaptic potentials (EPSPs) at p synapses and inhibitory postsynaptic potentials (IPSPs) at q inhibitory synapses. Once K(t) crosses Kthre from below, a spike is generated and K is reset to Krest. This model is termed as the IF model.LetKrest=0,a=b,p=q, and use the usual approximation to approximate the IF models (see Feng et al. [3] and Zhang and Feng [9]); then (1) can rewriten as(3)dKt=-1νKtdt+aλtdt+aλαtdWt,where {W(t)}t≥0 is a standard Brownian motion, α>0 is a constant, and if α=0.5, it implies that the input is derived from a Poisson process. If α>0.5, it is the so-called the supra-Poisson inputs, and the other is the so-called sub-Poisson inputs if α<0.5. In addition, a larger α leads to more randomness for the synaptic inputs.


## 2.2. General Linear Stochastic Differential Equation

In this subsection, we extend the one-dimensional I & F model (3) to n-dimensional stochastic differential equations in which the solution process enters linearly. Such processes arise in estimation and control of linear systems, in systems, in economics, and in various other fields (see Liu [10]), as(4)dXt=AtXt+Λtdt+ΣtdWt,X0=ξ,where W is an r-dimensional Brown motion independent of the n-dimensional initial vector ξ, and n×n, n×1, and n×r matrices A(t), Λ(t), and Σ(t) are nonrandom, measurable, and locally bounded, respectively.Now we define ann×n matrix function Φ(t), satisfying the following matrix differential equation:(5)Φ˙t=AtΦt,Φ0=I,where I is the d×d identity matrix. We know that (3) has unique (absolutely continuous) solution defined for 0≤t<∞, and, for each t≥0, the matrix Φ(t) is nonsingular.By Itô’s rule, it is easily verified that(6)Xt=Φtξ+∫0tΦ-1sΛsds+∫0tΦ-1sΣsdWs;0≤t<∞solves (4).We suppose thatEξ2<∞ and introduce the mean vector and covariance matrix functions as follows:(7)mt=EXt,ρs,t=EXt-msXt-mtT,Vt=ρt,t.From (4), we can show that(8)mt=Φtm0+∫0tΦ-1sΛsds,(9)ρs,t=ΦsV0+∫0s∧tΦ-1uΣuΣuΦ-1u′duΦ′t hold for every 0≤s,t<∞. In particular, m(t) and V(t) solve the linear equations:(10)m˙t=Atmt+Λt,V˙t=AtΛt+ΛtA′t+ΣtΣ′t.


## 3. Optimization Problem Formulation

We consider a simple model of high-order linear stochastic dynamical system ind-dimensional space. For simplicity of notation, we suppose that each component of trajectory in d-dimensional space satisfies the following stochastic differential equation (SDE):(11)xnt=anxn-1t+an-1xn-2t+⋯+a1xt+cλtdt+λtαdwt,where wt is a 1-dimensional Brown motion, and c>0,α>0,a1,…,an are constants. Generally we call λ(t) as the control signal. Now let(12)x1t=xtx2t=x˙t⋮xnt=xn-1t.Then we have (13)dx1t=x2tdtdx2t=x3tdt⋮dxnt=a1x1+a2x2+⋯+anxndt+cλtdt+λtαdwt.Thus, x1(t) is the position along some direction in space; the (n×n) matrix A(t) is a constant matrix:(14)A=010⋯0001⋯0⋮⋮⋮⋮⋮a1a2a3⋯anΛ(t)=c(0,0,…,λ(t))′ is the control signal vector, and Σ(t)=c0,0,…,λtα′ is a n×1 vector.Optimization Problem.For simplicity of notations, we let ξ=0. For a point dx∈R1 (in fact, dx expresses the arriving position component of the trajectory in the direction of component x(t)) and two positive numbers T,R, we intend to find a control signal λ∗(t) which satisfies the constrained condition:(15)x1t=dx,fort∈T,T+Rand such that the variance of x1(t) arrives the minimum in [T,T+R]; that is,(16)Iλ∗=minλ∈L2αT,T+RIλ=minλ∈L2αT,T+R∫TT+Rvarx1tdt.Let Φ(t)=(ϕij)n×n; by (6), we have, for 0≤t<∞,(17)x1t=c∫0tϕ1nt-sλsds+c∫0tϕ1nt-sλsαdws.Therefore, by (15) and (17), we have(18)∫0tϕ1nt-sλsds=dxc,t∈T,T+R.By the calculation of matrix Φ(t), we easily get the following results.Lemma 1.
Fork=1,2,…,n-1 and j=1,2,…,n,(19)ϕk+1,jt=ϕk,j′t,ϕk+1,jt=ϕ1,jkt.In particular, for n≥2,(20)ϕ1n0=0,ϕ1n′0=0,…,ϕ1nn-20=0.Proof.
SinceΦ′(t)=AΦ(t), by the definition of A (see (14)) and the multiplication of matrices, we get (19) at once.
SinceΦ(0)=I,Φk(0)=Ak, it is easy to get (20).
For simplicity of notation, we suppose thatA has n-different eigenvalues γ1,γ2,…,γn (in this case, A is diagonalizable). Hence A is similar to the diagonal matrix diag{γ1,γ2,…,γn} and there exists an invertible matrix S such that Φ(t)=S diag{exp(γ1t),…,exp(γnt)}S-1. Therefore, there are 1≤m≤n,m nonzero real number θ1,…,θm and m different eigenvalues γk1,…,γkm such that(21)ϕ1nt=∑i=1mθiexp⁡γkit≔∑i=1mθiexp⁡γit.We introduce the following notations:(22)Γt=diagexp⁡-γ1t,…,exp⁡-γmt;Θ=diagθ1,θ2,…,θm;Vγ1,…,γm=111⋯1γ1γ2γ3⋯γm⋮⋮⋮⋮⋮γ1m-1γ2m-1γ3m-1⋯γmm-1;Πt,λ=∫0texp⁡-γ1sλsds,…,∫0texp⁡-γmsλsds′D=dxc,0,…,0′.Theorem 2.
Under the constrained control condition (18), the following results hold:(i)
λ(t)=-a1dx/c,for t∈[T,T+R];(ii)
Π(T,λ)=Γ(T)Θ-1V-1(γ1,…,γm)D.Proof.
(i) SinceR>0, we can make the derivation until (n-1)-order derivation for the two sides of (17) and, by Lemma 1, we get for t∈[T,T+R](23)c∫0tϕ1nt-sλsds=dx∫0tϕknt-sλsds=0,k=2,…,n.That is to say,(24)∫0tΦt-sΛsds=cD,fort∈T,T+R.Since R>0,Φ(0)=I and Φ′(t)=AΦ(t), again differentiating (24), we obtain(25)Λt=-cAD∀t∈T,T+R. Therefore,(26)λt=-a1dxc,∀t∈T,T+R.
(ii) By (21) and Lemma 1, (23) can be expressed as(27)∑i=1mθiexp⁡γit∫0texp⁡-γisλsds=dxc,∑i=1mθiγik-1exp⁡γit∫0texp⁡-γisλsds=0,k=2,…,n;that is,(28)Vγ1,…,γmΘΓ-tΠt,λ=D,fort∈T,T+R. Since θ1,…,θm are different, using the multiplicity of matrix and T replacing t in the above equation, we get result (ii) at once.Note. if γ has multiplicity m>1 as an eigenvalue of A, and A is diagonal matrix, we can also choose n independent functions with the form tkexp⁡(γt),k=0,1,…,n-1. In this case, we can obtain the same result as that in the Theorem 2 by using the similar approach.By (17), it is easily seen that(29)Iλ=∫TT+Rvarx1tdt=c2∫TT+R∫0tϕ1n2t-sλs2αdsdt=c2∫0T∫TT+Rϕ1n2t-sdtλs2αds+c2∫TT+R∫sT+Rϕ1n2t-sdtλs2αds.Thus we only need to minimize the first term in (29), since minimizing each term of last equal in (29) implies minimizing I(λ) and, by (26), the control signal λs2α in the second term of (29) is a constant for s∈[T,T+R]. Now we apply the calculus of variations to the first term in (29), that is,(30)minλt∈L2α0,T+r∫0T∫TT+Rϕ1n2t-sdtλs2αds.To this end, let us define the control signal set: (31)λ·:∫0Tλsϕ1nT-sds=dxc,λt=-a1dxc,t∈T,T+R=UD,φ∈φ:∫0Tφsexp⁡-γisds=0,fori=1,2,…,m;φt=0,t∈T,T+R=UD0. For any τ∈R1,φ∈UD0 and λ∈UD, we have λ+τφ∈UD. By (30), λ∗ must satisfy (32)dIλ+τφdττ=0=0 which gives(33)∫0T∫TT+Rϕ1n2t-sdtλs2α-1signλsφsds=0.Comparing (33) with the first part constraint in UD0, we conclude that(34)∫TT+Rϕ1n2t-sdtλs2α-1signλs=∑i=1nξiexp⁡-γisalmost surely for s∈[0,T] with parameters ξi∈R1,i=1,2,…,n. Hence the solution of the original problem is(35)λ∗s=∑i=1nξiexp⁡-γis1/2α-1sign∑i=1nξiexp⁡-γis∫TT+Rϕ1n2t-sdt,fors∈0,Tλ∗s=-a1dxc,fors∈T,T+R,where ξ1,ξ2,…,ξn are the unique solution of the following (by the result of Theorem 2):(36)ΠT,λ∗=ΓTΘ-1V-1γ1,…,γmD.The similar equations are true for the other components ind-dimensional space.From the results above, we can obtain the following conclusions.Theorem 3.
Under the optimal control framework as we set up here andα>1/2, the optimal mean trajectory is a straight line. When α≤1/2 the optimal control problem is degenerate; that is, the optimal control signal is a delta function, and (29) with Λ=Λ∗ gives us an exact relationship between time T and variance.Proof.
This proof is similar to that of Theorem  1 in Feng et al. [3]; we omit it.Remark 4.
Whend=1, the results of Theorem 3 are consistent with Feng et al. [3]. The finding is also in agreement with the experimental Fitts law (see Fitts [11]); that is, the longer time of a reaching movement, the higher the accuracy of arriving at the target point.


## 4. Example in 3-Dimensional Space

We consider a simple model of (arm) movements. LetX(t)=(x(t),y(t),z(t))′ be the position of the hand at time t; we then have(37)X¨=-1τ1τ2X-τ1+τ2τ1τ2X˙+1τ1τ2Λtdt+Λα,tdW→t,where τ1, τ2, α>0 are parameters, Λ(t)=(λx(t),λy(t),λz(t))′ is the control signal, Λ(α,t) is a 3×3 diagonal matrix with diagonal elements as λztα,λytα,λxtα, respectively, and W→(t)=(Wx(t),Wy(t),Wz(t))′ is the standard Brown motion. In physics, we know that (37) is the well-known Kramers’ equation. In neuroscience, it is observed in all in vivo experiments that the noise strength is proportional to the signal strength Λ(t) and hence the signals received by muscle take the form of (37) (see Feng et al. [3]).For a pointD=(dx,dy,dz)′∈R3 and two positive numbers T, R, we intend to find a control signal Λ(t) which satisfies(38)μt=Xt=D,fort∈T,T+R,(39)IΛ∗=minλ∈L2αT,T+RIΛ=minλ∈L2αT,T+R∫TT+REXt-μt′Xt-μtdt=minλ∈L2αT,T+R∫TT+Rvarxt+varyt+varztdt,where Λ∈L2[0,T+R] means that each component of it is in L2[0,T+R]. To stabilize the hand, we further require that the hand will stay at D for a while, that is, in time interval [T,T+R], which also naturally requires that the velocity should be zero at the end of movement. The physical meaning of the problem we considered here is clear; at time T, the hand will reach the position D (see (38)), as precisely as possible (see (39)). Without loss of generality, we assume that dx>0,dy>0, and dz>0.To use the results in the previous section, we can rewrite the optimal control problem posed in the previous paragraph as a2-order linear stochastic dynamical system in 2-dimensional space, that is,(40)x¨t=-1τ1τ2xt-τ1+τ2τ1τ2x˙t+1τ1τ2λxtdt+λxtαdWxt.The similar equation holds true for y(t) and z(t).If we letvx(t) express the moving velocity in the direction of x-coordinate, (40) becomes the following 2-order linear SDE:(41)dxt=vxtdt,dvxt=-1τ1τ2xt-τ1+τ2τ1τ2vxt+1τ1τ2λxtdt+1τ1τ2λxtαdWxt.Comparing (12), it is easy to know that (42)A=01-1τ1τ2-τ1+τ2τ1τ2,Λt=01τ1τ2λt,Σt=01τ1τ2λxtα.Since(43)A=S-1τ100-1τ2S-1,where (44)S=-τ1-τ211,S-1=1τ2-τ11τ2-1-τ1,thus by calculating, we know that(45)etA=ϕ11tϕ12tϕ21tϕ22t=1τ2-τ1τ2exp⁡-tτ2-τ1exp⁡-tτ1τ1τ2exp⁡-tτ2-exp⁡-tτ1exp⁡-tτ1-exp⁡-tτ2τ2exp⁡-tτ1-τ1exp⁡-tτ2,where(46)ϕ12t=τ1τ2τ2-τ1exp⁡-tτ2-exp⁡-tτ1,ϕ22t=τ1τ2τ2-τ11τ1exp⁡-tτ1-1τ2exp⁡-tτ2.Hence, by (8), we get that(47)xt=1τ1τ2∫0tϕ12t-sλsds+1τ1τ2∫0tϕ12t-sλsαdBxs,vxt=1τ1τ2∫0tϕ22t-sλsds+1τ1τ2∫0tϕ22t-sλsαdBxs.Therefore, (34), (35), and (36) become (48)∫TT+Rϕ122t-sdtλs2α-1signλs=ξ1exp⁡sτ1+ξ2exp⁡sτ2,almost surely for s∈[0,T] with parameters ξi∈R,i=1,2,…,n. Hence the solution of the original problem is(49)λ∗s=ξ1exp⁡s/τ1+ξ2exp⁡s/τ21/2α-1signξ1exp⁡s/τ1+ξ2exp⁡s/τ2∫TT+Rϕ1n2t-sdt,with ξ1,ξ2 being given by the following equations (by Theorem 2):(50)∫0Tλ∗sexp⁡sτ1ds=τ1dxexp⁡Tτ1,∫0Tλ∗sexp⁡sτ2ds=τ2dxexp⁡Tτ2.


## 5. Conclusion

The experimental study of movement in human has shown that voluntary reaching movements obey Fitts law: the longer the time taken for a reaching movement, the greater the accuracy for the hand to arrive at the end point. In this paper, we study a stochastic control problem for a reaching movement within ad-dimensional space. We solve this stochastic control problem explicitly and obtain the analytical solutions for optimal signals, optimal velocity, and optimal variance. Furthermore, we find that the optimal control is also consistent with Fitts law. This implies that the straight line trajectory is a natural consequence of optimal stochastic control principles, under a nondegenerate optimal control signal.



---

*Source: 8730859-2017-12-17.xml*