Electrical Impedance Tomography (EIT) is a non-invasive medical imaging technique in which an image of the electrical properties (conductivity or permittivity) of part of the body is inferred from surface electrode measurements. It has the advantages over other techniques like X-rays and requires no exposure to radioactive materials. Applications include detection of breast cancer, pulmonary emboli, blood clots, impaired gastric emptying and etc.. Essentially, through only surface measurements, the internal electric conductivity and permittivity are identified as an image inside the human body. For instance, the electric conductivity of malignant tumor, a high-water-content tissue, is one order higher than that of the normal (fat) tissue, which allows one to find potential diseases and locations through the constructed image throughout the body .
EIT is also a useful tool in other fields such as geophysics, environmental sciences and nondestructive testing of materials. It is able to locate underground mineral deposits, detect leaks in underground storage tanks and monitor flows of injected fluids into the earth for extraction or environmental cleaning. Moreover, EIT can detect the the corrosion or defects of construction material and machine parts   when invasion testing is not possible or destructive.
For the applications mentioned above, researchers are often faced with the problem of how to work out the conductivity inside an object when only part of the boundary measurements are available . It is well known that this inverse problem is nonlinear, unstable and intrinsically ill-posed .
In theory, complete boundary measurements indeed determine the conductivity in the interior uniquely  , however, in practice only limited number of electrodes and current patterns are available from measurements. Various numerical algorithms to reconstruct the conductivity have been proposed and fall into two categories, non-iterative and interative methods. Noniterative methods were developed based on the assumption that the conductivity does not differ too much from a constant. Calderón  proved that a map between the conductivity constant and a quadratic energy functional in (4) is injective when is close to a constant in a sufficiently small neighborhood, and an approximation formula was also given to reconstruct the conductivity. The back-projection method of Barber-Brown  gave a crude approximation to the conductivity increment based on inverse of the generalized Radon transform, which works best for smooth or whose singularity is far from the boundary. For a L-electrode system, Noser algorithm minimizes the sum of squares error of the voltages on the electrodes by using one step of a regularized Newton’s method. Meanwhile, the iterative methods are devoted to minimize different regularized least squares functionals such as a Tikhonov-type regularization  and a total variation formulation  where iterative gradient-based optimization algorithms are commonly used.
Solving inverse problems with iterative algorithms usually requires a solution of a forward problem at each iteration numerically, the computation time accumulates fast for commonly used grid based global methods such as FEM/BEM methods. In most of the EIT problem, the measurements are only available on limited number of electrodes, in finding the conductivity profile to match the measured voltages there, a global solution of the potentials over the whole object is in fact not needed, therefore a global solution procedure during hte forward problem incur unneccessary computational cost beyond the electrodes. With this in mind, in this paper, an alternative local stochastic approach based on path integral Monte Carlo (PIMC) simulations using Feynman-Kac formula will be proposed. Due to the nature of the Feynma-Kac formula which allows the potential solution at any single location including those on the electrodes, we could dramatically reduce the amount of solutions needed for each forward problem solution.
The remainder of the paper is organized as follows. In Section 2, the forward and inverse problems of the EIT problem are introduced. Section 3 includes the path integral Monte Carlo method for the forward mixed boundary value problem of the voltage-to-current map for a 8-electrode EIT problem for a spherical object. Also, in order to validate the accuracy of the stochastic method, we include a boundary element method to generate reference solutions. Comparison between the stochastic and deterministic methods shows the accuracy of the proposed method. Finally conclusions are drawn in Section 4.
2 Forward and inverse problems in EIT
2.1 The forward problem
In this section, we will first review the forward problem arising from EIT. The mathematic models for EIT have been developed and compared with the experimental measurement of voltages on electrodes for a given conductivity distribution, which is adjusted to fit the measurements. The existing models are continuum model, gap model, shunt model and complete electrode model. Among all, the complete electrode model was shown to be capable of predicting the experimentally measured voltages to within 0.1 percent  and the existence and uniqueness of the model has also been proved .
Complete Electrode Model
Let the domain of the object is denoted as embedded within we assume there is an anomaly . The domain is assumed to have a smooth boundary with a limited number of electrodes attached to . The conductivity inside is given by and the electric potential for the model will satisfy the following boundary value problem equationparentequation
|In (1f), we have prescribed a constant potential on the surface of the tumor, which corresponds to modeling the highly conductive anomaly as a perfect conductor.|
Equation (1a) describes the distribution of electric potential in the interior of the object where is the conductivity or inverse of the resistivity. It was derived from Maxwell equations by neglecting the time-dependence of the alternating current and assuming the current source inside the object to be zero . Both equations (1b) and (1c) indicates the knowledge of current density on and off the electrodes on the boundary, respectively. Equation (1d) takes account of the electrochemical effect by introducing
as the contact impedance or surface impedance which quantitatively characterizes a thin, highly resistive layer at the contact between the electrode and the skin, which causes potential jumps according to the Ohms law. It should be noted that, the regularity of potentialdecreases as the contact impedance approaches zero , which becomes a huge hindrance to accurate numerical resolution as in practice usually good contacts with small contact impedance are used.
The usual choice of numerical method to solve the forward problem is grid based methods such as finite element method, boundary element method. The global solution will require the solution of a large linear system as fine meshes ususally are needed around the electrodes. Our approach in this paper is to develop a local stochastic method, i.e. a path integral Monte Carlo method based on on Feynman-Kac formula for the mixed boundary value problem in (1a)-(1d). This will allows us to only find the unknown solution on a single electrode only without solving any global matrix system, thus dramatically reducing the cost of a forward problem solver in an iterative inverse problem solution of EIT.
2.2 The inverse problem
The EIT inverse problem proposed by Calderón in 1980 posed the question if it was possible to determine the heat conductivity of static temperature and heat flux measurements on the boundary? The inverse problem can be formulated mathematically as follows.
Let be a bounded domain in , with a Lipschizian boundary , and be a real bounded measurable function in with a positive lower bound. Consider the differential operator
acting on and a quadratic form for functions in
The problem is then to decide whether is uniquely determined by the operator and to calculate in terms if is indeed determined by . To put it in another way, we need to verify the injectivity of the following map
In the context of physical meaning, represents the power necessary to maintain an electrical potential on . Calderoń showed that the map is analytic if and is injective. He also proved that determines when is sufficiently close to a constant.
Literally speaking, can be determined through boundary measurements as we will see below. Then the problem is reduced to whether the conductivity can be reconstructed through the surface electrode measurements. Many authors have made contributions to the problem under various assumptions. Kohn and Vogelius  provided a positive answer to the determination of the conductivity which is in and has all derivatives at the boundary. Sylvester and Uhlmann  proved uniqueness for conductivities in the plane while Brown  relaxed the regularity of the conductivity to derivatives.
We may treat Calderón Problem from another perspective. To be specific, knowing for each is equivalent to the knowledge of “Dirichlet-to-Neumann” data. In fact, by the Green’s identity and assuming that the potential on the boundary of the
If setting above, then we have
where the left-hand side of (6) is exactly and the right-hand side involves Neumann values if Dirichlet conditions are given on the boundary. Thus, the Calderón Problem can be restated as whether is uniquely determined by the “Dirichlet-to-Neumann” map on the boundary. The map tells us how the boundary potential determines the current flux across the boundary . Also it is clear that a “Robin-to-Neumann” map (or voltage- to-current map) is essentially equivalent to a “Dirichlet-to-Neumann” map through a close look at (1d). Given a full Robin-to-Neumann map , it uniquely determines and thus equivalent to the Dirichlet-to-Neuman map. Under such circumstances, uniqueness of solutions to the inverse conductivity problem was proved by Astala and Päivärinta  without any regularity imposed on the boundary for a bounded measurable conductivity in two dimensions. In three dimensions, Haberman and Tararu  confirmed the answer for conductivities and Lipchitz conductivities close to the identity.
The Robin-to-Neumann map is also called a voltage-to-current map. Under the Complete Electrode model, to obtain the voltage-to-current map is equivalent to solving the forward problem on the whole boundary given a known conductivity. As discussed before, the traditional finite element method or/and boundary element method may require dense mesh near the contacts of skin and electrodes, resulting in large linear systems to be solved. Therefore, more efficient schemes are desired without constraints from the geometries of the domain and shapes of the electrodes.
One possible way to improve the efficiency would be to develop a probabilistic estimator of the voltage-to-current map. The main advantage of the method lies in the prevailing multicore computing. Maire and Simon proposed a so-called partially reflecting random walk on spheres algorithm to compute voltage-to-current map in a parallel manner, which is also efficient when only solutions at only a few points are desired. For the Dirichlet boundary problem, it is well known that killed Brownian motion is the stochastic process that governs the differential operator. However, in the mixed boundary situations, the partially reflecting Brownian motion comes into play in preventing the path running out of the domain by either absorption or instantaneous reflection. Simulation of absorption is much easier to take care of comparing to that of reflection as the latter requires special techniques like local finite difference discretization. Various schemes of first order or second order schemes have been proposed and analyzed 
. Maire and Simon also studied a similar approach involving second order space discretization scheme. A variance reduction technique was introduced as well to improve the efficiency and accuracy of the method.
In our work, we aim to find a probabilistic solution to the voltage-to-current map by directly simulating the reflecting Brownian motion paths on the boundary. More precisely, the calculation of the boundary local time is treated explicitly in details and then integrated into the Feynman-Kac type representation, and the voltages can be then obtained numerically on the boundary.
3 Numerical scheme for forward problem and voltage-to-current map
In medical applications, limited number of electrodes are attached to human body to get surface measurements. We will illustrate our numerical method for a model problem for a unit spherical object. In Fig. 4, eight electrodes are superimposed on the boundary and the centers of the electrodes all lie on the plane with radius 0.2.
Consider the conductivity equation (conductivity is taken to be 1 outside the anomaly without loss of generality) with both Neumann and Robin boundary value problem
where is a constant between 0 and 1 and is the outward unit normal verctor.
In our previous work , we described a method based on Monte Carlo simulations to find potentials on boundaries through WOS sampling. The same approach is employed to find voltages on the electrode patches given mixed boundary conditions for the Laplace operator. From the Robin boundary conditions on the electrode patches, the Neumann data are automatically known, thus, the voltage-to-current map is obtained for the forward EIT problem.
3.1 A path integral Monte Carlo Solution using Feynman-Kac formula
First, let us review some preliminaries concerning Feynman-Kac formula for the mixed boundary value problem in (7), which lays the foundation of our path integral Monte Carlo approach.
3.1.1 Reflecting Brownian Motion and boundary local time
Assume that is a domain with a boundary in . A generalized
Skorohod problem is stated as follows:
Let , a continuous function from to . A pair is a solution to the Skorohod equation if
is continuous in ;
the local time is a nondecreasing function which increases only when , namely,
The Skorohod equation holds:
denotes the outward unit normal vector at.
The Skorohod problem was first studied in  by A.V. Skorohod in addressing the construction of paths for diffusion processes inside domains with boundaries, which experience the instantaneous reflection at the boundaries. Skorohod presented the result in one dimension in the form of an Ito integral and Hsu  later extended the concept to -dimensions ().
In general, solvability of the Skorohod problem is closely related to the smoothness of the domain . For higher dimensions, the existence of (10) is guranteed for domains while uniqueness can be achieved for a domain by assuming the convexity for the domain . Later, it was shown by Lions and Sznitman  that the constraints on can be relaxed to some locally convex properties.
Suppose that is a standard Brownian motion (SBM) path starting at and is the solution to the Skorohod problem , then will be the standard reflecting Brownian motion (SRBM) on starting at
. Because the transition probability density of the SRBM satisfies the same parabolic differential equation as that for a BM, a sample path of the SRBM can be simulated simply as that of the BM within the domain. However, the zero Neumann boundary condition for the density of SRBM implies that the path be pushed back at the boundary along the inward normal direction whenever it attempts to cross the boundary.
The boundary local time is not an independent process but associated with SRBM and defined by
where is a strip region of width containing and . Here is the local time of , a notion invented by P. Lévy . This limit exists both in and -. for any .
It is obvious that measures the amount of time that the standard reflecting Brownian motion spends in a vanishing neighborhood of the boundary within the time period . An interesting part of (11) is that the set has a zero Lebesgue measure while the sojourn time of the set is nontrivial . This concept is not just a mathematical one but also has physical relevance in understanding the “crossover exponent” associated with “renewal rate” in modern renewal theory .
3.1.2 Simulation of RBM and calculation of local time
Method of WOS for Brownian paths
Random walk on spheres (WOS) method was first proposed by Müller , which can solve the Dirichlet problem for the Laplace operator efficiently.
To illustrate the WOS method for the Dirichlet problem of Laplace equation, with Dirichlet boundary conditions . The solution can be rewritten in terms of a measure defined on the boundary ,
where is the harmonic measure defined by
It can be shown easily that the harmonic measure is related to the Green’s function for the domain with a homogeneous boundary condition , i.e.,
If the starting point of a Brownian motion is at the center of a ball, the probability of the BM exiting a portion of the boundary of the ball will be proportional to the portion’s area. Therefore, sampling a Brownian path by drawing balls within the domain can significantly reduce the path sampling time. To be specific, given a starting point inside the domain , we simply draw a ball of largest possible radius fully contained in
and then the next location of the Brownian path on the surface of the ball can be sampled, using a uniform distribution on the sphere, say at. Treat as the new starting point, draw a second ball fully contained in , make a jump from to on the surface of the second ball as before. Repeat this procedure until the path hits a absorption -shell of the domain (see Fig. 2) . When this happens, we assume that the path has hit the boundary (see Fig. 1(a) for an illustration).
Now we can define an estimator of (12) by
where is the number of Brownian paths sampled and is the first hitting point of each path on the boundary. To speed up the WOS process, maximum possible size of the sphere for each step would allow faster first hitting on the boundary, see Fig. 1(b).
For the reflecting boundary, we will construct a strip region around the boundary (see Fig. 2) and allow the process to move according to the law of BM continuously. Before the path enters the strip region, the radius of WOS is chosen to be of a maximum possible size within the distance to the boundary. Once the particle is in the strip region, the radius of the WOS sphere is then fixed at a constant (or , see Fig. 3). With this approach, according to the definition (11), the local time may be interpreted as
given a prefixed constant in the strip region and be the cumulative steps that path stays within the -region from the begining until time (see Remark below for definition). Notice that only those steps where the path of remains in the -region will contribute to while the SRBM may lie out of the -region at other steps. More details can be found in . One may refer to Fig. 3 for an illustration of the behavior of path near the boundary.
Remark Occupation time of SRBM in the numerator of was calculated in terms of that of BM sampled by the walks on spheres. Notice here that within the -region, the radius of the WOS may be or , which implies that the corresponding elapsed time of one step for local time could be or . The latter is four times bigger than the former. But if we absorb the factor into , still holds. In practical implementation, we treat as a vector of entries of increasing value, the increment of each component of over the previous one after each step of WOS will be 0, 1 or 4, corresponding to the scenarios that is out of the -region, in the -region while sampled on the sphere of a radius , or in the -region while sampled on the sphere of a radius , respectively.
3.1.3 Feynamn-Kac formula
We consider the mixed boundary value problem in the domain to the mixed problem
The probabilistic solution for the boundary value problem above is given by the well-known Feynman-Kac formula 
where is the standard reflecting Brownian motion, is the corresponding local time and the Feynman-Kac functional , is the first time a Brownian path originating from hits the boundary of .
Feynman-Kac formula provides a local solution procedure to solve the partial differential equations through stochastic processes. As a matter of fact, the infinitesimal generator of the Laplace operator is Brownian motion which is involved in the solution as well. For a general elliptical operator, Itó processes come into play.
The numerical approximation to will be
where denotes each step of the path and denotes the steps where the path hits the boundary (Robin or Neumann).
In the context of complete electrode model, the second expectation is zero due to the zero Neumann boundary. One may find more details in .
3.2 A deterministic solution with the boundary element method
To provide reference solutions to the path integral MC method proposed above, we will present a deterministic method based on boundary element method for the mixed boundary value problem (7).
3.2.1 Graded boundary mesh
Let be the unit ball centered at the origin and , the collection of electrodes is . Each electrode patch is assumed to have an equal surface area.
Boundary mesh is constructed on the electrode patches and off-electrodes, respectively as shown in Fig.5 and Fig.6. For our implementation, GMSH is used to generate an unstructured 2D mesh consisting of flat triangles given a “size field” while on the electrode patches the meshes are structured in such a way that mesh points are found at the intersection of division along the longitude and altitude which gives a body-fitted mesh and the mapping from the elemental triangle to the curved ones can be found in . A global boundary integral equation can then be set up based on the two meshes.
Since the contact impedance varies from 0 to 1, we take for our numerical tests. A close look at the boundary conditions in (7) reveals that discontinuities at the rims of all the electrodes. It is natural to enlarge the radius of mesh on so that we may have an easy control over the mesh size for calculation. Assume the enlarged radius to be . Because of the discontinuity, we consider a graded mesh on the enlarged surface by introducing a layered mesh structure, as Fig. 7 illustrates. There are four layers: the first ranges from center to , second from to , third from to and fourth from to . A dense mesh will be used around the rim () of the electrode, which implies that both 2nd layer and 3rd layer should have a decreased mesh size towards . Furthermore, a graded mesh also discretizes the first layer while an evenly distributed mesh is used on the fourth layer. And and are the number of divisions along the altitude in each layer, respectively. Here we take and . The mesh size can be calculated through . The number of divisions along longitude will be the same for each layer, i.e. . Fig. 8 shows the realization of the graded mesh on the north pole patch. The red points are the mesh points on the electrode (first two layers) and blue ones are off-electrode points on the rest two layers. We can see clearly the mesh on can be constructed similarly or obtained through rotation of that on along -axis. Besides, the “size field” of GMSH is 0.012 on D which yields 88383 mesh points and 175530 flat triangle elements as Fig. 6 shows.
3.2.2 Boundary integral equation
By using the second Green’s identity, we have for any fixed ,
where and is a Green’s function for the exterior domain of with homogeuous boundary condition on , i.e.,
As from the interior, using the fact that both and vanish on the boundary we obtain equationparentequation
where the Robin condition in (7) is used to deduce from (27a) to (27b). Therefore we obtain a global boundary integral equation for . With the different meshes on and , we let sweep over all the mesh points where (28) is imposed, we have a linear system for solving on . As a result, reference potentials on the boundary can be obtained to validate the Monte Carlo simulations. Meanwhile, any reference potentials inside the domain can be achieved as well through
where , the Neumann values on the boundary, are automatically known once the reference potentials are found on the electrode patches from the Robin conditions.
3.2.3 Reference current on electrodes
With the preparations in the previous sections, the results are shown in terms of current on each electrode according to
A direct summation of over eight electrodes yields the whole current to be -7.10e-7 which is close to 0 as the conservation of charges condition suggests. Meanwhile the electrode currents show symmetric patterns with respect to both and axis, consistent with the system design.
3.3 Validation of path integral MC numerical results
All Monte Carlo simulation is based on the algorithm in (24). The number of Monte Carlo simulations is for all the mesh points on the electrodes while the length of path ranges from 1900 to 3000 on different patches. From Table 2, the absolute errors between the numerical approximations and reference currents remain at a low level not higher than 1% .
Walk-on-sphere method is used to simulate the reflecting Brownian motion within the spherical object. Feynman-Kac formula shows an infinite length of path while from the numerical perspective it is truncated to be a finite number in such a way that is increased until the desired accuracy is achieved and it is lowered when is increased further through our test.
With the numerical approximations of potentials on the boundary, the Neumann values are automatically known and thus a full Robin-to-Neumann map has been achieved over the whole boundary by conducting the same procedure at all the mesh points on different electrode patches. Therefore, we’ve found an effective way to compute the voltage-to-current map without resorting to the conventional deterministic methods which involves constructing 3-D mesh in the domain and solution of a large global matrix. It has also advantages when the geometry of the domain in 3-D is complicated.
4 Conclusions and future work
This paper presents a path integral Monte Carlo method to solve the forward problem of EIT and the voltage-to-current map is acquired as needed for the iterative algorithm of the inverse EIT problem where the forward problem needs to be solved first at each iteration. The method takes advantage of the parallel computing capability in modern multicore computers and solutions can be calculated simultaneously at different electrodes independently.
In our calculations, the contact impedance is taken to be 0.5 while it can be very close to zero. In that case, the voltages will jump drastically and then a dense mesh must be placed at the contacts otherwise the calculation of the reference voltages may be undermined. However, the PIMC will not be affected as each point is independently calculated.
Meanwhile, local boundary integral equations may be considered as an alternative approach to the voltage-to-current map. To be more specific, when the reference potentials are obtained through the global boundary integral equation, then potentials at any interior points are known due to (29). Then, we may construct local boundary integral equations of Neumann values on each electrode.
The illness of the inverse problem requires high accuracy of forward modelling otherwise it may lead to very large fluctuations in reconstructions of conductivity. Our method is accurate to some extent but needs more development in convergence and variance reduction.
W.C. acknowledges the support of the National Science Foundation (DMS-1764187) for the work in this paper.
-  Borcea, Liliana, Electrical impedance tomography, Inverse problems 18.6 (2002): R99.
-  Calderón , Alberto P., On an inverse boundary value problem, Computational & Applied Mathematics 25.2-3 (2006): 133-138.
-  Astala, Kari, and Lassi Pivrinta, Caldern’s inverse conductivity problem in the plane, Annals of Mathematics (2006): 265-299.
-  J. A. Given, Chi-Ok Hwang and M. Mascagni, First-and last-passage Monte Carlo algorithms for the charge density distribution on a conducting surface, Physical Review E 66, 056704, 2002.
-  J. P. Morillon, Numerical solutions of linear mixed boundary value problems using stochastic representations, Int. J. Numer. Meth. Engng., Vol. 40, 387-405, 1997.
-  A.V. Skorokhod, Stochastic equations for diffusion processes in a bounded region, Theory of Probability & Its Applications 6.3 (1961), 264-274.
-  H. Tanaka, Stochastic differential equations with reflecting boundary condition in convex regions, Hiroshima Mathematical Journal 9.1 (1979), 163-177.
-  P.L. Lions and A.S. Sznitman, Stochastic differential equations with reflecting boundary conditions, Communications on Pure and Applied Mathematics 37.4 (1984), 511-537.
-  I. Karatzas and S. E. Shreve, Brownian motion and stochastic calculus, Springer-Verlag New York Inc., 1988.
-  J.F. Douglas, Integral equation approach to condensed matter relaxation, Journal of Physics: Condensed Matter 11.10A (1999), A329.
-  M. E. Müller, Some continuous Monte Carlo methods for the Dirichlet problem, The Annals of Mathematical Statistics, Vol. 27, No. 3, 569-589, 1956.
-  Borsic, Andrea, Brad M. Graham, Andy Adler, and William RB Lionheart, Total variation regularization in electrical impedance tomography, (2007).
-  P. Lévy, Processus Stochastiques et Mouvement Brownien, Gauthier-Villars, Paris, 1948.
-  Alessandrini, Giovanni, Stable determination of conductivity by boundary measurements, Applicable Analysis 27.1-3 (1988): 153-172.
-  Santosa F, Vogelius M, A backprojection algorithm for electrical impedance imaging, SIAM Journal on Applied Mathematics, (1990)216-43.
-  R. Kohn and M. Vogelius, Determining conductivity by boundary measurements, Comm. Pure Appl. Math., 37 (1984), pp. 113-123.
-  J. Sylvester and G. Uhlmann, A uniqueness theorem for an inverse boundary value problem in electrical prospection, Comm. Pure Appl. Math., 39 (1986), pp. 91-112
-  Kaipio, Jari P., et al, Statistical inversion and Monte Carlo sampling methods in electrical impedance tomography, Inverse problems 16.5 (2000): 1487.
-  Somersalo, Erkki, Margaret Cheney, and David Isaacson, Existence and uniqueness for electrode models for electric current computed tomography, SIAM Journal on Applied Mathematics 52.4 (1992): 1023-1040.
-  Vauhkonen, Marko, et al., Tikhonov regularization and prior information in electrical impedance tomography, IEEE transactions on medical imaging 17.2 (1998): 285-293.
-  Maire, Sylvain, and Martin Simon, A partially reflecting random walk on spheres algorithm for electrical impedance tomography, Journal of Computational Physics 303 (2015): 413-430.
-  Russell M. Brown, Global Uniqueness in the Impedance-Imaging Problem for Less Reg- ular Conductivities, SIAM Journal on Mathematical Analysis 27 (1996), no. 4, 1049.
-  Haberman, Boaz, and Daniel Tataru, Uniqueness in Calderan’s problem with Lipschitz conductivities, Duke Mathematical Journal 162.3 (2013): 497-516.
-  Audus, Debra J., et al, Interplay of particle shape and suspension properties: a study of cube-like particles, Soft matter 11.17 (2015): 3360-3366.
-  Cheney, Margaret, David Isaacson, and Jonathan C. Newell., Electrical impedance tomography, SIAM review 41, no. 1 (1999): 85-101.
-  Juba, Derek, Walid Keyrouz, Michael Mascagni, and Mary Brady, Acceleration and Parallelization of ZENO/Walk-on-Spheres, Procedia Computer Science 80 (2016): 269-278.
-  S. Maire and E. Tanré, Monte Carlo approximations of the Neumann problem, Monte Carlo Methods and Applications 19.3 (2013): 201-236.
-  Bossy, Mireille, et al, Probabilistic interpretation and random walk on spheres algorithms for the Poisson-Boltzmann equation in Molecular Dynamics ESAIM: Mathematical Modelling and Numerical Analysis 44.5 (2010): 997-1048.
-  Lin, Huimin, Huazhong Tang, and Wei Cai, Accuracy and efficiency in computing electrostatic potential for an ion channel model in layered dielectric/electrolyte media, Journal of Computational Physics 259 (2014): 488-512.
-  K. L. Chung, Green, Brown, and Probability and Brownian Motion on the Line, World Scientific Pub Co Inc, 2002.
-  (Elton) P. Hsu, “Reflecting Brownian motion, boundary local time and the Neumann problem”, Dissertation Abstracts International Part B: Science and Engineering[DISS. ABST. INT. PT. B- SCI. ENG.], Vol. 45, No. 6, 1984.
-  Y. Zhou, W. Cai and (Elton) P. Hsu, Computation of Local Time of Reflecting Brownian Motion and Probabilistic Representation of the Neumann Problem, Commun. Math. Sci., Vol. 15(2017), 237-259.
-  Y. Zhou, and W. Cai., Numerical Solution of the Robin Problem of Laplace Equations with a Feynman-Kac Formula and Reflecting Brownian Motions, Journal of Scientific Computing 69.1 (2016): 107-121.