Multiscale Modeling of Light Absorption in Tissues: Limitations of Classical Homogenization Approach

In biophotonics, the light absorption in a tissue is usually modeled by the Helmholtz equation with two constant parameters, the scattering coefficient and the absorption coefficient. This classic approximation of “haemoglobin diluted everywhere” (constant absorption coefficient) corresponds to the classical homogenization approach. The paper discusses the limitations of this approach. The scattering coefficient is supposed to be constant (equal to one) while the absorption coefficient is equal to zero everywhere except for a periodic set of thin parallel strips simulating the blood vessels, where it is a large parameter The problem contains two other parameters which are small: , the ratio of the distance between the axes of vessels to the characteristic macroscopic size, and , the ratio of the thickness of thin vessels and the period. We construct asymptotic expansion in two cases: and and prove that in the first case the classical homogenization (averaging) of the differential equation is true while in the second case it is wrong. This result may be applied in the biomedical optics, for instance, in the modeling of the skin and cosmetics.


Physical background of the problem
In the present paper we consider the Helmholtz equation with rapidly oscillating large potential with a periodic support having a small measure of its intersection with a period. This absorption coefficient q of the potential (more precisely, it is the ratio between the absorption coefficient and the diffusion coefficient) depends on three small parameters: E is a standard homogenization parameter that is the ratio of the period of the potential and the characteristic macroscopic size; d is the ratio between the measure of the intersection of the support of q and the period; v {1 is a small parameter standing for the inverse of the ratio of the maximal value of the coefficient q and the diffusion coefficient multiplied by the square of the characteristic macroscopic size of the problem (we will consider the case when q takes only two values: v and 0).
The Helmholtz equation is considered below as a model of the light absorption in tissues under hypothesis that this absorption takes place only in the set of parallel thin blood vessels (where q~v=0) and this absorption is ignored outside of these vessels. The linear dimensions of the vessels are much smaller than the linear dimensions of the body as a whole. Typically the distance between two neighboring large micro-vessels (of about 15 mm in the diameter) is around 180 mm in primate cerebral cortex [1].
This distance is also used in 3D simulation of tumor growth and angiogenesis [2]. Assume at the first approximation that the tissue is a nearly periodic structure. Moreover, consider the twodimensional idealization of this periodic structure that is, the periodic set of the parallel strait narrow identic vessels separated by the homogeneous tissue. Let L be the macroscopic characteristic size and let EL be the distance between two neighboring vessels (strips). It means that the parameter E stands here for the ratio of the distance between two neighboring vessels and the characteristic macroscopic size L. It is assumed throughout that this ratio is a small parameter. Indeed, if the characteristic macroscopic size (L) is equal to 10 mm (it is a typical value of the diffuse optical tomography [3]) and the distance between the vessels is equal to 0.18 mm, then E~0:018. The second parameter of the model is the ratio of the thickness of the vessels and the distance between neighboring vessels denoted d. This d as well is supposed to be a small parameter. Thus, in the discussed above structure we have d~15=180~0:08333 . . .. In the spectral window 500 nm-700 nm, the oxyhaemoglobin extinction coefficient shows a wide dynamic from 275 (mole/L) {1 cm {1 at 690 nm (the lowest) and 55 500 (mole/L) {1 cm {1 at 575 nm [4]. In order to quantify the maximum variation we compare arterial blood (90 percent saturated) with pure water at 690 nm and 575 nm then the ratio of the absorption coefficient is respectively 300 and 400000. At 700 nm the absorption of tissue (without blood) is about ten times less than that of the pure water [5]. In the visible window the maximum of this ratio can be more than 10000. That is why in the idealized model we consider the case when the absorption coefficient is equal to zero out of vessels. As we have mentioned above the light absorption process is described here by the Helmholtz equation (1), where q is just nondimensionalized absorption coefficient, equal to zero out of vessels and equal to the great (dimensionless) parameter v within the vessels. Let us remind that the physical sense of this great parameter is the ratio of the absorption and diffusion effects, more exactly, the ratio of the absorption coefficient and the diffusion coefficient multiplied (i.e. the ratio is multiplied) by the square of the characteristic macroscopic size L. Let the haemoglobin concentration be 150 g/liter. Then in order to convert the molar extinction coefficient E to the absorption coefficient (in mm {1 ), one has to multiply it by 0.00054. At the wavelength 575 nm for the scattering coefficient close to 2.5 mm {1 , v is about 23000. Then E 2 vd~0:62 and E 2 vd 2~0 :05: In the case of penetrating vessels with the average diameter close to 65 mm, at the same wavelength, E 2 vd 2~1 . Then the case when E 2 vd 2 is greater than one can be found in a spectral window below 600 nm and for vessels of large diameter. In neurophotonics, the highly vascularized pie-matter correspond to the maximal values of E 2 vd 2 . The highly vascularized tumors could be also a special case when E 2 vd 2 is very high.
Models of this ''composite medium'' are widely considered in biophysics. Let us give a short review of these results. The partial differential equation (PDE) diffusion with absorption is one of the most important in life science, and in particular, in tissue optics [6]. The researchers are interested in two types of results: a) when the optical properties are unknown, and they want to find them using observed measurements, it is the harder wellknown inverse problem. b) when the optical properties are known, then they can use many strategies to calculate various quantities of interest (reflectance, transmission, fluence). This is often referred to as the forward problem.
But in all these approaches, the community of research in biophotonics uses the volumetric averaging of the absorption coefficient. This problem is crucial in the domain of in vivo diffuse optical tomography [3,7], and is always present in in vivo optical neuromethods [8].
There are as well some experimental papers on the problem of averaging of blood absorption of light in tissue but without the homogenization theory [9][10][11][12][13][14].
Our goal in the present paper is to show that the classical homogenization (averaging) approach to the solution of equation (1) leading to the approximation where vqw is often approximated to the volumic mean value of q, has some limitations. Indeed, we will show that it is right for some combination of magnitudes of parameters E,d It means that the ''diluting'' of vessels in multiscale modeling of light absorption of tissues may be applied only in justified cases, in particular, in case A. In fact, the situations when the coefficients of equations depend on two or more small parameters and it leads to a fail of classical homogenization are not too surprising now: first results of this type could be found in [15][16][17]. For some cases the non-classical asymptotic expansions were constructed. However, to our knowledge, such examples have not been constructed for rapidly oscillating large potentials with narrow support. The present paper tries to contribute in this important case of biophysical applications.
In order to simplify technical details of the analysis we consider the boundary value problem for the Helmholtz equation set in a layer 0vx 1 v1,x 2 [R f gwith the Dirichlet conditions on the boundary of the layer; we assume that the right-hand sides of the equation and of the boundary condition do not depend on x 2 . In this case we may seek a solution independent of x 2 as well. The Helmholtz equation takes the following form: where x stands for x 1 , E is a small positive parameter, such that, 1=E is an integer. We consider a boundary condition corresponding to a constant solution in the case of the constant absorption coefficient q~v and of a constant right-hand side f~1. Then u E,d,v~1 =v. So, these boundary conditions are: Similar problems were addressed by many authors [18][19][20][21][22]; however the potentials considered in these works are typically E {1 or E {2 . In this work, there are three parameters and we study all relevant asymptotic regimes. We follow the ideas and methods of [15][16][17]23,24] in the sense that we construct asymptotic expansions for u E and analyze all the important asymptotic regimes of the parameters E, v, d.
Mention that the Helmholtz equation (1) (or (3)) is not a perfect model for the light absorption process: it is not more than an approximation of a more adequate model based on the radiation transfer equation [25]. Indeed, one can introduce the average (over all directions v) for the diffused radiation intensity I(x,v) : In the case of the described above plane geometry the diffusion approximation of the radiation transfer has a form [25] { d dx 1 where [(x 1 ) and s(x 1 ) are the coefficients corresponding to the light absorption and dispersion respectively, g is the blood heat radiation intensity and f is the external radiation intensity, L is the thickness of the irradiated skin strip, 1 3 is the dimension factor.
Introducing the new variable corresponding to the relative optical depth of the radiation process we get for U(X ) (here U(X (x 1 ))~u(x 1 )): Making the change of unknown function we get finally equation (3): The choice of the boundary conditions (4) is not too important for our analysis: the asymptotic approach developed below can be applied in the case of other conditions, for instance, corresponding to a constant solution for equation (3) with q replaced by its volumic mean vqw, or or the periodicity conditions (mention that the diffusion approximation of the radiation transfer may loose its precision near the boundary).

Mathematical statement of the problem
Consider interval V~½0,1 and let 0vdv1=2: Let q(:) be a 1periodic function defined on the basic period ½0,1 by Let f be smooth enough (this assertion will be formulated more precisely later). We are interested in studying the asymptotic behaviour as E?0, v??, and d?0 of u E,d,v solving the following boundary value problem: Note that we have the following a priori estimate for solution of equation (6) with boundary condition (7) replaced by the homogeneous one Proposition 0.1. There exists a constant C independent of small parameters such that Here C may be taken equal to 1=2. A well-known imbedding theorem gives: Corollary 0.1. There exists a constant C independent of small parameters such that Applying equation, we get: Corollary 0.2. There exists a constant C independent of small parameters such that Returning now to equation (6) with non-homogeneous boundary condition (7) we can present the solution as a sum of the solution of equation (6) with homogeneous boundary condition and the solution of the homogeneous equation (6) with nonhomogeneous boundary condition (7). Applying the maximum principle argument, we get: Proposition 0.2. There exists a constant C independent of small parameters such that solution of problem (6)-(7) satisfies: Remark 0.1. If boundary condition (7) is replaced by another one: with some numbers A and B then the estimate of Proposition 0.2 takes form: Let us mention one more useful inequality for functions v of Proof: For any x, If v 0 (x 0 )~0, then the estimate is proved. If v 0 (x 0 )=0, then the sign of v 0 is constant on the whole interval ½0,l; consider the case v 0 w0: Then v is increasing and v(1)~v(0)z Ð l 0 v 0 (t)dt §v(0)zlv 0 (x 0 ), and so, and we get (8). In the same way we consider the opposite case when v 0 v0:

Methods
Asymptotic expansion in the case E?0,v??,d?0,vd??,E 2 vd?0 Consider the K-th level approximation u K E of u E given by an ansatz from [15]: where N p l (:) are 1-periodic functions, D l~d l dx l , and v E is a smooth function which will be sought in a form Substituting the expression for u K E from (9) in (6) we obtain where H p l~d 2 N p l dj 2 (j 1 )z2 and B E,K 1 is given by In writing expressions (11), (12) we followed the convention that N p l~0 if at least one of the two indices p or l, is negative.
We choose H p l (j) being a constant such that the equation (12) has a solution, i.e.
where S.T~ð .dj: Also, periodic solutions N p l of equation are chosen to be of average zero for all p and l such that, (p,l)=(0,0); N 0 0 is chosen equal to 1. We get: N 0 0~1 , N 0 l~0 , N 1 l~0 for all l, N 2 0 is a 1-periodic solution of equation

SqT~vd:
We solve the family of equations (15)  Proof. The proof is by induction on m.
Step (i). Let m~0. That is, p~0,1. We get: N 0 0~1 , h p l~S N p l{2 (j)T~d l,2 d p,0 , N p l satisfy equations (15) with the last two terms equal to zero; so N 0 l~0 , for all lw0, N 1 l~0 for all l §0 and the assertion of lemma for m~0 is evident.
Step (ii). Let m §1 and let all N p l , D j N p l and h p l with pƒ2mz1 be bounded by c l,m (vd) ½p=2 , where ½a stands for the integer part of a, and c l,m is a constant independent of the small parameters. Consider equation (15) for p~2mz2. For l~0 it takes form with the piecewise polynomial right hand side; moreover, the primitive (with vanishing mean value) of this right hand side is also bounded by the same order. That is why, integrating this equation twice and keeping every time the vanishing mean for the primitive, we prove the assertion of lemma for l~0. In the same way we prove it for N 2mz3 0 . Now we apply again (15), express D 2 j N p l from this equation and prove by induction on l that all h 2mz2 l , h 2mz3 l , the right hand sides and so, finally, N 2mz2 l and N 2mz3 l are bounded by C l,mz1 (vd) mz1 . The lemma is proved.
Substituting the expression for u K E from (9) in (7) we obtain ð13Þ and applying the 1{periodicity of functions N p l and the divisibility of 1 by E we get: Substituting the expression for v E (10) in (11), (16),(17), we get: and Here and Define now functions v r from boundary value problems: and Here SqT~dv is supposed to be greater than 1. These equations have a form where a r ,b r are some real numbers and f r is a smooth function. In the same way as in the first section we prove that where the constants C and C 1 are independent of small parameters.
Applying the induction on m, we get: where the constants C m are independent of small parameters. Applying now (8), we get: for mw0 and so, Remark 0.2. Functions v r depend on small parameters but their asymptotic expansion can be constructed by classical boundary layer technique (see [23]). For instance, if f [C 2K1 (½0,1), then v 0 has a form: Moreover, if f [C 2K1zK2 (½0,1), then R 0 is K 2 times differentiable and satisfies the estimate (see the above a priori estimate for v r ): In the same way we construct by induction the expansions of functions v r , rƒKz1. Assume that K 1 §2K, K 2 §Kz1. Consider the right hand side f r having the form of a linear combination with some bounded coefficients h p l .
(vd) ½p=2 of functions g having a form: where F 0r,j ,F 1r,j ,F r,j are 2K 1 zK 2 {r times differentiable functions independent of small parameters vd and E, such that, there exist constants c 1 ,c 2 independent of small parameters satisfying for any real j inequalities D l j F 0r,j (j) ƒc 1 e {c 2 j j j , D l j F 0r,j (j) ƒc 1 e {c 2 j j j , lƒK 2 {l, F r,j are K 1 zK 2 {r times differentiable on ½0,1, The right hand sides of the boundary conditions as well have a similar form: they are some linear combinations with bounded coefficients N p l (0) .
(vd) ½p=2 of constants s and t having a form: where a r,j ,b r,j are independent of small parameters.
The expansion of v r has a similar form of a linear combination with bounded coefficients of functions v rg (x)X where u r,j satisfy equations (functions with the negative indices are equal to zero) and exponentially decaying at infinity functions U 0r,j satisfy equations .
Applying the induction on r, the explicit expressions for the right hand sides of the problems (24), (25), (26) for v r and the estimates of Lemma 0.1 we prove that there exist constants C l,r independent of small parameters, such that, This estimate and Lemma 0.1 are crucial for evaluation of the discrepancies B 1 ,B 2 ,B 3 and B 4 : Now, applying the a priori estimate of the Remark after Proposition 0.2 we get the estimate Assume that there exists a positive real c such that, v(E 2 vd) c is bounded by a constant independent of small parameters. Then the last bound yields: Taking K 1 w2czK, we get: On the other hand, (27) and Lemma 0.1 give: So, from the triangle inequality we get Theorem 0.1. Assume that there exists a positive real c such that, v(E 2 vd) c is bounded by a constant independent of small parameters. Let f belong to C 3(2czKz1) (½0,1). Then there exists a constant C independent of small parameters such that Consider now the very first term of the expansion v 0 and K~0 we get: Corollary 0.3. There exists a constant C independent of small parameters such that Mention that the error of order E 2 is much smaller than v 0 that is of order 1 vd , and so the relative error of the approximation of the exact solution u E,d,v by v 0 is small. Remark 0.3. Here we have considered the case when vd stands for a large parameter tending to the infinity. The case when vd is a positive finite constant may be considered by a classical technique [15] and the estimate of Corollary 0.3 becomes Asymptotic analysis: case E?0, v??, d?0, E 2 vd 2 ??
Here we will construct an example where the behavior of the solution is completely different from the behavior described in the previous section. In particular, the solution of the problem is different from the behavior of the solution v 0 of the homogenized equation. For the sake of simplicity consider the right hand side f (x)~1. Consider an auxiliary function g defined on the interval ½{1,1 as a function from C 2 (½{1,1) independent of small parameters such that it equals to zero on the interval ½{1=3,1=3 and it equals to one on ½{1,{2=3|½2=3,1. Let us keep the same notation g for the 2-periodic extension of this function on R. Consider an approximation for u E,d,v having a form: theoretically. The homogenized model may be inapplicable in the case if E 2 vd is not small.
On the effective absorption coefficient, the volumetric mean and the total absorption coefficient Biological tissues are highly heterogeneous media at the microscopic scale. For years, the patterns of the mammalian cortex micro-angioarchitecture have attracted the interest of many research groups [1,[26][27][28]. At different scales the blood vessels have very different physical characteristics. Moreover the absorption coefficient is great in the vessels and small out of the vessels in UV and visible parts of the spectrum. At the end they induce multiscale complex photon propagation and absorption. On the other hand the experimental data are obtained mainly at the macroscopic scale, and so one of the main questions is how to make a passage from the micro-scale to the macro-scale. The main mathematical tool of such passage from one scale to another (the up-scaling) is the homogenization theory, see [15,19,23,24] and the references there. Normally, the microscopic description of the heterogeneous medium can be replaced by an equation with constant coefficients. This equation is called the homogenized equation and the constant coefficients are effective coefficients. This homogeneous approximation for the heterogeneous medium is justified if the solutions of both models are close in some norms. In some cases a heterogeneous medium cannot be approximated by a homogeneous one [16,17,23,24]. Then the notion of an effective coefficient cannot be introduced in the above sense.
In particular, in the present paper we prove that if product E 2 vd is small then the homogenized model (2) is justified and the effective absorption coefficient is equal to the volumetric mean value vqw of function q. The smallness of the value of this product E 2 vd means that the relative error of the approximation of the heterogeneous medium by the homogeneous one is of order of this product. The effective absorption coefficient is an important quantity because the detailed knowledge of the tissue macroscopic optical properties is essential for an optimization of optical methods i.e. for modeling the color of skin, of port-wine stains and for tumor detection; it helps to adapt an appropriate photodynamic therapy, in particular, some laser treatment.
On the other hand the present paper shows the limitations of the homogenized models for the absorption problems. Indeed, if product E 2 vd 2 is not small then the homogenized model (2) as well as other homogenized models (with other possible constant values of the effective absorption coefficient) do not approximate the initial microscale model, because the solution of problem (2) does not oscillate for any choice of the effective absorption coefficient, while the solution of equation (6) with boundary condition (7) rapidly oscillates (see (29)).
This theoretical argument is confirmed by some physical reasons and by experimental observations. Many authors [10,12,14,29] show that the assumption of a homogeneous distribution of blood in the tissue may strongly overestimate the total blood absorption when absorption is high and/or the vessels have sufficient diameter. For large vessels less of light reaches the center of the vessel, and the absorbers in the center of a vessel contribute less and less to the total attenuation of the light. However, the total blood absorption is an important qualitative characteristic of the absorption process and so it should be calculated with a great precision. Let us apply the above asymptotic analysis provided for equation (6) with boundary condition (7) and calculate the total blood absorption in cases A and B. Let us define the total absorption coefficient m as a ratio of the integrals For the particular case of a constant coefficient q we get evidently m~q. In the case A substituting the leading term of an asymptotic expansion we get that the asymptotic behavior of m is given by an approximate formula m&vqw, i.e. it is close to the volumetric mean of q that is, vd. In the case B the value of this total absorption coefficient m is much less than the volumetric mean (that was observed in the discussed above experiments): the direct computations for the approximation (29) show that m&vqw=(0:5E ffiffiffi ffi v p d) with E ffiffiffi ffi v p dww1: The authors of papers [10,12,14,29] characterize the fall of the total absorption coefficient by the correction factor C d~m =vqw and discuss the situations when this factor is different from 1. In our asymptotic analysis we see that in the case A C d &1 but in the case B C d vv1.
We hope that our result will help to analyze the link between the total absorption coefficient m and other parameters which may be applied in optical tomography [3,7].