CN104637085A - Building method of decoupling fluorescence monte-carlo model - Google Patents

Abstract

The invention relates to a building method of a decoupling fluorescence monte-carlo model. The method comprises the following steps that an emitting density integral equation is built and is used for describing the transmission of exciting light in biological tissues and the fluorescence exciting and transmission process; through the core function conversion in the fluorescence exciting and fluorescence transmission process, optical parameters of the fluorescence exciting and transmission process are decoupled to be correlated with the path probability density function, so that a sampling fluorescence photon path is identical to a laser photon path, and meanwhile, a weight function of fluorescence photons is correlated with the optical parameters in the fluorescence exciting and transmission process; the weight of the fluorescence photons is directly calculated according to the weight function of the fluorescence photons. The model built by the method has the characteristics that the precision is high, and the model can be applicable to any turbid medium and any fluorescence distribution, and the method has wide application prospects in the aspects of fluorescence imaging, fluorescence spectrum measurement and analysis.

Description

A kind of method for building up separating coupling fluorescence monte-Carlo model

Technical field

The present invention relates to mathematical simulation and biomedical engineering field, be specifically related to a kind of method for building up separating coupling fluorescence monte-Carlo model.

Background technology

Imaging-PAM and fluorescence spectral measuring occupy critical role in biomedical sector, and it is widely used in cancer diagnosis, medicament research and development and the visual research of gene expression.Along with fluorescent probe technique significant progress, Imaging-PAM start to be applied to the large molecule activity rule of the inner specific biological of small animal model in volume tracing and measurement.Under the effect of outside excitation source, Internal Fluorescent group of biological tissue absorbs the light of certain wavelength and sends fluorescence, and the fluorescence spectrum detected can reflect the function information of molecule.In recent years, fluorescent molecular tomography technology obtains and pays close attention to widely.

In order to improve the counting yield of fluorescent molecular tomography, many research groups all focus on fluorescence Monte Carlo method and accelerate.Liebert proposes and applies in stratiform turbid media, efficient time-resolved fluorescence Monte Carlo method [1].In this approach, by preserving a large amount of photon path information, the analytic relationship of transmission photons weight and tissue optical parameter is established.Jin Chen proposes a kind of high-efficiency fluorescence Monte Carlo method [2] based on gate fluorescence Jacobian matrix computing time.Background weight matrix from excitaton source to detector can be drawn by photon path and fluorescent absorption coefficient calculations.Anand T.N.Kumar proposes a kind of absorption perturbation Monte Carlo method based on preserving photon path information, is applicable to the fluorescent sample [3] of time resolution multiparameter.Above three kinds of methods, in fluorescence Monte Carlo simulation process, have all done some supposition as identical by the path of anisotropic scattering replacement exciting light and fluorescence random walk in the tissue in isotropic scatterning during fluorescence excitation.In order to fluorescent molecular tomography is quantitative better, it is very crucial for setting up that a high precision and high efficiency computation model go to describe photon transmission in biological tissues.Above three kinds of methods have done more supposition in simulation process, and therefore the usable condition of method receives a definite limitation, therefore cannot meet the demand of fluorescent molecule tomography rebuilding completely.Therefore invent a kind of by the Kernel Function Transformation in fluorescence excitation and fluorescence communication process, solution is coupled out fluorescence excitation and the optical parametric of communication process and associating of path probability density function, make the fluorescent photon path of sampling consistent with excitation light subpath, make the weighting function of fluorescent photon be associated with the optical parametric of fluorescence excitation and communication process simultaneously; Thus along excitation light subpath, directly calculate the fluorescence monte-Carlo model of the weight of fluorescent photon according to fluorescence weighting function, be applicable to the situation of any fluorescence distribution of any turbid media, and meet this demand.

[1]Liebert A，Wabnitz H，Obrig H，et al.Non-invasive detection of fluorescence from exogenous chromophores in the adult human brain[J].Neuroimage，2006，31(2)：600-608.

[2]Chen J，Venugopal V，Intes X.Monte Carlo based method for fluorescence tomographic imaging with lifetime multiplexing using time gates[J].Biomedical optics express，2011，2(4)：871-886.

[3]Kumar A T N.Direct Monte Carlo computation of time-resolved fluorescence in heterogeneous turbid media[J].Optics letters，2012，37(22)：4783-4785.

Summary of the invention

The object of the invention is by the Kernel Function Transformation in fluorescence excitation and fluorescence communication process, solution is coupled out fluorescence excitation and the optical parametric of communication process and associating of path probability density function, make the fluorescent photon path of taking out consistent with excitation light subpath, make the weighting function of fluorescent photon be associated with the optical parametric of fluorescence excitation and communication process simultaneously; Thus along excitation light subpath, directly calculate the weight of fluorescent photon according to fluorescence weighting function.The model computational accuracy that the present invention sets up is high and be applicable to the situation of any fluorescence distribution of any turbid media.

Separate a method for building up for coupling fluorescence monte-Carlo model, it is characterized in that comprising the following steps:

Step 1: build Particle transport equations and describe photon transmitting procedure in biological tissues under stable state;

Step 2: the emission density integral equation building transport equation equivalence describes process that exciting light photon transmits in biological tissues and exciting light photon is absorbed and inspires fluorescent photon, and fluorescent photon continues the process transmitted in biological tissues;

Step 3: the kernel function in fluorescence excitation process in the neumann series solution of integral equation and fluorescence communication process is all represented by the kernel function of exciting light communication process, thus solution is coupled out fluorescence excitation and the optical parametric of communication process and associating of path probability density function;

Step 4: throw in exciting light photon, follows the trail of photon transmission in biological tissues, utilizes exciting light scattering coefficient to sample, along exciting light path, and the fluorescence intensity on calculating detector.

In step 1, under stable state, photon transmitting procedure is in biological tissues described as:

$\begin{array}{c}\hat{s}\·\▿\mathrm{\φ}(\stackrel{\→}{r},\hat{s},v,t)+{\mathrm{\μ}}_t(\stackrel{\→}{r},v,t)\mathrm{\φ}(\stackrel{\→}{r},\hat{s},v,t)=S(\stackrel{\→}{r},\hat{s},v,t)+\\ {\∫}_{4\mathrm{\π}}{\mathrm{\μ}}_t(\stackrel{\→}{r},v,t)\mathrm{\φ}(\stackrel{\→}{r},{\hat{s}}^{\′},v,t)\·C({\hat{s}}^{\′}\→\hat{s}){\mathrm{d\Ω}}^{\′}\end{array}$

for the direction of photon current state;

for the direction of the next state of photon current state;

R is the position of photon current state;

V is the frequency of photon current state;

The solid angle unit of the next state that d Ω ' is photon current state;

for in t, near r point, in cell frequency interval and unit solid angle unit, direction exists neighbouring photon number;

for being photon in t, be emitted in the neighbouring cell frequency interval of r point and unit solid angle unit by external source, direction exists neighbouring photon number;

for direction is photon, at r point, place enters collision, and collision rift direction exists the photon probability number in direction;

μ _tfor extinction coefficient.

The process that in step 2, exciting light transmission and fluorescence excitation also transmit is described as respectively:

$x\left(p\right)=\∫x\left(p^{\′}\right)K(p^{\′}\→p,{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}}+{\mathrm{\μ}}_{\mathrm{af}})d{p}^{\′}+S\left(p\right)$

$y\left(p\right)=\∫y\left(p^{\′}\right)K(p^{\′}\→p,{\mathrm{\μ}}_s^{\mathrm{em}},{\mathrm{\μ}}_a^{\mathrm{em}}){\mathrm{dp}}^{\′}+\∫x\left(p^{\′}\right)K\left(p^{\′}\right)K(p^{\′}\→p,{\mathrm{\μ}}_{\mathrm{af}}){\mathrm{dp}}^{\′}$

P is photon transmission state, and p ' is the photon transmission state of after p.

The emission density that the emission density that x (p) is exciting light, y (p) are fluorescence.

S (p) is source emission density.

K is the kernel function of light communication process.

for exciting light scattering coefficient, for the absorption coefficient of exciting light.

for fluorescence scattering coefficient, for the absorption coefficient of fluorescence.

μ _affor the absorption coefficient of fluorophore.

In step 3, solution is coupled out associating of fluorescent optics parameter (exciting and communication process) and path probability density function, and the kernel function of all available exciting light communication process of kernel function in fluorescence excitation process and fluorescence communication process represents:

(1) kernel function of the kernel function exciting light communication process of fluorescence excitation process is expressed as:

$\begin{array}{c}K(p^{\′}\→p,{\mathrm{\μ}}_{\mathrm{af}})=T({\stackrel{\→}{r}}^{\′}\→\stackrel{\→}{r}|{\hat{s}}^{\′},{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}}+{\mathrm{\μ}}_{\mathrm{af}})C({\hat{s}}^{\′}\→\hat{s}|\hat{r},{\mathrm{\μ}}_{\mathrm{af}})\\ =K({\stackrel{\→}{r}}^{\′}\→\stackrel{\→}{r}|{\hat{s}}^{\′},{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})\exp (-{\∫}_0^{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}{\mathrm{\μ}}_{\mathrm{af}}({\stackrel{\→}{r}}^{\′}+l\hat{s})\mathrm{dl})\×\\ \frac{{\mathrm{\μ}}_t^{\mathrm{ex}}\left(\stackrel{\→}{r}\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(\stackrel{\→}{r}\right)}{{\mathrm{\μ}}_t^{\mathrm{ex}}\left(\stackrel{\→}{r}\right)}\frac{C({\hat{s}}^{\′}\→\hat{s}|\stackrel{\→}{r},{\mathrm{\μ}}_{\mathrm{af}})}{C({\hat{s}}^{\′}\→\hat{s}|\stackrel{\→}{r},{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})}\\ =K(p^{\′}\→p,{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})\exp (-{{\∫}_0^{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}\mathrm{\μ}}_{\mathrm{af}}({\stackrel{\→}{r}}^{\′}+l\hat{s})\mathrm{dl})\×\\ \frac{\mathrm{\η}{\mathrm{\μ}}_{\mathrm{af}}\left(\stackrel{\→}{r}\right)}{{\mathrm{\μ}}_s^{\mathrm{ex}}\left(\stackrel{\→}{r}\right)}\frac{P_I({\hat{s}}^{\′}\·\hat{s})}{P_A({\hat{s}}^{\′}\·\hat{s})}\end{array}$

The position of the next state that r ' is photon current state;

for the photon sent from r ', enter the density of collision at r place.

for the extinction coefficient of exciting light.

L is the step-length that photon often walks.

η is quantum efficiency.

P _ifor isotropic phase bit function, P _afor anisotropy phase function.

(2) kernel function of the kernel function exciting light communication process of fluorescence communication process is expressed as:

$\begin{array}{c}K(p^{\′}\→p,{\mathrm{\μ}}_s^{\mathrm{em}},{\mathrm{\μ}}_a^{\mathrm{em}})=T({\stackrel{\→}{r}}^{\′}\→\stackrel{\→}{r}|{\hat{s}}^{\′},{\mathrm{\μ}}_s^{\mathrm{em}},{\mathrm{\μ}}_a^{\mathrm{em}})C({\hat{s}}^{\′}\→\hat{s}|\stackrel{\→}{r},{\mathrm{\μ}}_s^{\mathrm{em}},{\mathrm{\μ}}_a^{\mathrm{em}})\\ =T({\stackrel{\→}{r}}^{\′}\→\stackrel{\→}{r}|{\hat{s}}^{\′},{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})\exp (-{\∫}_0^{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}({\mathrm{\μ}}_t^{\mathrm{em}}({\stackrel{\→}{r}}^{\′}+l\hat{s})-{\mathrm{\μ}}_t^{\mathrm{ex}}({\stackrel{\→}{r}}^{\′}+l\hat{s}))\mathrm{dl})\×\\ C({\hat{s}}^{\′}\→\hat{s}|\stackrel{\→}{r},{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})\frac{{\mathrm{\μ}}_s^{\mathrm{em}}\left(\stackrel{\→}{r}\right)}{{\mathrm{\μ}}_s^{\mathrm{ex}}\left(\stackrel{\→}{r}\right)}\\ =\frac{{\mathrm{\μ}}_s^{\mathrm{em}}\left(\stackrel{\→}{r}\right)}{{\mathrm{\μ}}_s^{\mathrm{ex}}\left(\stackrel{\→}{r}\right)}\exp (-{\∫}_0^{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}({\mathrm{\μ}}_t^{\mathrm{em}}({\stackrel{\→}{r}}^{\′}+l\hat{s})-{\mathrm{\μ}}_t^{\mathrm{ex}}({\stackrel{\→}{r}}^{\′}+l\hat{s}))\mathrm{dl})K(p^{\′}\→p,{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})\end{array}$

for the extinction coefficient of fluorescence.

In step 4, the fluorescent photon path of sampling is consistent with excitation light subpath, and the weighting function of fluorescent photon is associated with fluorescent optics parameter (exciting and communication process); Along excitation light subpath, directly calculate the weight of fluorescent photon according to fluorescence weighting function.

Fluorescence intensity in step 4 on detector is expressed as:

$\begin{array}{c}D={\∫y}_m\left(p\right)g\left(p\right)\mathrm{dp}\\ =\underset{m=0}{\overset{\∞}{\mathrm{\Σ}}}\∫...\∫{\mathrm{\τ}}_m\left(p\right)w_m\left(p\right)\underset{k=0}{\overset{m}{\mathrm{\Π}}}{\mathrm{dp}}_k\mathrm{dp}\end{array}$

D is the fluorescent photon intensity that detection receives.

The number of times that m collides for photon before arrival detector.

The function of fluorescence statistic of g (p) for detector detects.

τ _mp probability density that () is photon.

W _mp weighting function that () is photon.

(1) the probability density τ of photon _mp () is expressed as:

${\mathrm{\τ}}_m\left(p\right)=S\left(p_0\right)K(p_0\→p_1,{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})...K(p_{m-1}\→p_m,{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})$

Wherein S (p)=S (p ₀), S (p ₀) be source emission density;

(2) the weighting function w of photon _mp () is expressed as:

$\begin{array}{c}{w}_m\left(p\right)=g\left(p\right)\underset{i=0}{\overset{m}{\mathrm{\Σ}}}\exp (-{\∫}_0^{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}{\mathrm{\μ}}_{\mathrm{af}}({\stackrel{\→}{r}}^{\′}+l\hat{s})\mathrm{dl})\frac{{\mathrm{\η\μ}}_{\mathrm{af}}\left({\stackrel{\→}{r}}_i\right)}{{\mathrm{\μ}}_s^{\mathrm{ex}}\left({\stackrel{\→}{r}}_i\right)}\\ \×\frac{P_I({\hat{s}}_i\·{\hat{s}}_{i+1})}{P_A({\hat{s}}_i\·{\hat{s}}_{i+1})}\underset{j=i+1}{\overset{m}{\mathrm{\Π}}}\left(\frac{{\mathrm{\μ}}_s^{\mathrm{em}}\left({\stackrel{\→}{r}}_j\right)}{{\mathrm{\μ}}_s^{\mathrm{ex}}\left({\stackrel{\→}{r}}_j\right)}\right)\exp (-{\∫}_0^{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}({\mathrm{\μ}}_t^{\mathrm{em}}({\stackrel{\→}{r}}^{\′}+l\hat{s})-{\mathrm{\μ}}_t^{\mathrm{ex}}({\stackrel{\→}{r}}^{\′}+l\hat{s}))\mathrm{dl})\end{array}$

Accompanying drawing explanation

Fig. 1 is basic flow sheet of the present invention.

Fig. 2 is the cylinder model figure of concrete structure of the present invention.

Fig. 3 a is normalization fluorescence intensity profile on the detector that obtains in the simulation of standard fluorescence Monte Carlo (sfMC) method.

Fig. 3 b to separate in the simulation of coupling fluorescence Monte Carlo (dfMC) method normalization fluorescence intensity profile on the detector that obtains.

Fig. 4 is the level line comparison diagram that standard conciliates normalization light intensity distributions image on the detector that obtains in coupling fluorescence monte carlo modelling.

Embodiment

The invention will be further described by reference to the accompanying drawings.

As shown in Figure 1, implementation step of the present invention is as follows:

(1) build Particle transport equations and describe photon transmitting procedure in biological tissues under stable state;

$\begin{array}{c}\hat{s}\·\▿\mathrm{\φ}(\stackrel{\→}{r},\hat{s},v,t)+{\mathrm{\μ}}_t(\stackrel{\→}{r},v,t)\mathrm{\φ}(\stackrel{\→}{r},\hat{s},v,t)=S(\stackrel{\→}{r},\hat{s},v,t)+\\ {\∫}_{4\mathrm{\π}}{\mathrm{\μ}}_t(\stackrel{\→}{r},v,t)\mathrm{\φ}(\stackrel{\→}{r},{\hat{s}}^{\′},v,t)\·C({\hat{s}}^{\′}\→\hat{s}){\mathrm{d\Ω}}^{\′}\end{array}$

(2) the emission density integral equation building transport equation equivalence describes exciting light photon and transmits in biological tissues, exciting light photon is absorbed and is inspired fluorescent photon, and fluorescent photon continues to transmit the process being finally detected device and receiving in biological tissues;

Exciting light photon transmits emission density integral equation of equal value in biological tissues and is described as:

$x\left(p\right)=\∫x\left(p^{\′}\right)K(p^{\′}\→p,{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}}+{\mathrm{\μ}}_{\mathrm{af}})d{p}^{\′}+S\left(p\right)$

Exciting light photon is absorbed and is inspired fluorescent photon, and the process prescription that fluorescent photon continues to transmit in biological tissues is:

$y\left(p\right)=\∫y\left(p^{\′}\right)K(p^{\′}\→p,{\mathrm{\μ}}_s^{\mathrm{em}},{\mathrm{\μ}}_a^{\mathrm{em}}){\mathrm{dp}}^{\′}+\∫x\left(p^{\′}\right)K\left(p^{\′}\right)K(p^{\′}\→p,{\mathrm{\μ}}_{\mathrm{af}}){\mathrm{dp}}^{\′}$

(3) kernel function in fluorescence excitation process in the neumann series solution of integral equation and fluorescence communication process is all represented by the kernel function of exciting light communication process, set up contacting of fluorescence excitation process and fluorescence communication process and exciting light communication process, thus solution is coupled out fluorescence excitation and the optical parametric of communication process and associating of path probability density function;

Kernel function K represents the conversion from state p to p ', can be expressed as the product of collision kernel and migration core:

$K(p^{\′}\→p,{\mathrm{\μ}}_s,{\mathrm{\μ}}_a)=T({\stackrel{\→}{r}}^{\′}\→\stackrel{\→}{r}|{\hat{s}}^{\′},{\mathrm{\μ}}_s,{\mathrm{\μ}}_a)C({\hat{s}}^{\′}\→\hat{s}|\stackrel{\→}{r},{\mathrm{\μ}}_s,{\mathrm{\μ}}_a)$

T is migration core, represents the photon sent from r ', enters the density of collision at r place.

$T({\stackrel{\→}{r}}^{\′}\→\stackrel{\→}{r}|{\hat{s}}^{\′},{\mathrm{\μ}}_s,{\mathrm{\μ}}_a)={\mathrm{\μ}}_t\left(\stackrel{\→}{r}\right)\exp (-{\∫\mathrm{\μ}}_t({\stackrel{\→}{r}}^{\′}+l\hat{s})\mathrm{dl})\mathrm{\δ}(\hat{s}-(\stackrel{\→}{r}-+{\stackrel{\→}{r}}^{\′})/|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}\left|\right)/{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}^2$

C is collision kernel, represents that direction is photon, at r point, place enters collision, and collision rift direction exists photon probability number in the d Ω solid angle unit in direction.

The collision kernel of exciting light communication process is expressed as:

$C({\hat{s}}^{\′}\→\hat{s}|\stackrel{\→}{r},{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})={\mathrm{\μ}}_s^{\mathrm{ex}}\left(\stackrel{\→}{r}\right)P_A({\hat{s}}^{\′}\·\hat{s})/({\mathrm{\μ}}_t^{\mathrm{ex}}\left(\stackrel{\→}{r}\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(\stackrel{\→}{r}\right))$

The collision kernel of fluorescence excitation process is expressed as:

$C({\hat{s}}^{\′}\→\hat{s}|\stackrel{\→}{r},{\mathrm{\μ}}_{\mathrm{af}})={\mathrm{\μ}}_{\mathrm{af}}\left(\stackrel{\→}{r}\right)P_A({\hat{s}}^{\′}\·\hat{s})/({\mathrm{\μ}}_t^{\mathrm{ex}}\left(\stackrel{\→}{r}\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(\stackrel{\→}{r}\right))$

The collision kernel of fluorescence communication process is expressed as:

$C({\hat{s}}^{\′}\→\hat{s}|\stackrel{\→}{r},{\mathrm{\μ}}_s^{\mathrm{em}},{\mathrm{\μ}}_a^{\mathrm{em}})={\mathrm{\μ}}_s^{\mathrm{em}}\left(\stackrel{\→}{r}\right)P_A({\hat{s}}^{\′}\·\hat{s})/({\mathrm{\μ}}_t^{\mathrm{em}}\left(\stackrel{\→}{r}\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(\stackrel{\→}{r}\right))$

In non-fluorescence region-type be 0.

The kernel function of the kernel function exciting light communication process of fluorescence excitation process is expressed as:

$\begin{array}{c}K(p^{\′}\→p,{\mathrm{\μ}}_{\mathrm{af}})=T({\stackrel{\→}{r}}^{\′}\→\stackrel{\→}{r}|{\hat{s}}^{\′},{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}}+{\mathrm{\μ}}_{\mathrm{af}})C({\hat{s}}^{\′}\→\hat{s}|\hat{r},{\mathrm{\μ}}_{\mathrm{af}})\\ =K({\stackrel{\→}{r}}^{\′}\→\stackrel{\→}{r}|{\hat{s}}^{\′},{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})\exp (-{\∫}_0^{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}{\mathrm{\μ}}_{\mathrm{af}}({\stackrel{\→}{r}}^{\′}+l\hat{s})\mathrm{dl})\×\\ \frac{{\mathrm{\μ}}_t^{\mathrm{ex}}\left(\stackrel{\→}{r}\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(\stackrel{\→}{r}\right)}{{\mathrm{\μ}}_t^{\mathrm{ex}}\left(\stackrel{\→}{r}\right)}\frac{C({\hat{s}}^{\′}\→\hat{s}|\stackrel{\→}{r},{\mathrm{\μ}}_{\mathrm{af}})}{C({\hat{s}}^{\′}\→\hat{s}|\stackrel{\→}{r},{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})}\\ =K(p^{\′}\→p,{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})\exp (-{{\∫}_0^{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}\mathrm{\μ}}_{\mathrm{af}}({\stackrel{\→}{r}}^{\′}+l\hat{s})\mathrm{dl})\×\\ \frac{\mathrm{\η}{\mathrm{\μ}}_{\mathrm{af}}\left(\stackrel{\→}{r}\right)}{{\mathrm{\μ}}_s^{\mathrm{ex}}\left(\stackrel{\→}{r}\right)}\frac{P_I({\hat{s}}^{\′}\·\hat{s})}{P_A({\hat{s}}^{\′}\·\hat{s})}\end{array}$

The kernel function of the kernel function exciting light communication process of fluorescence communication process is expressed as:

$\begin{array}{c}K(p^{\′}\→p,{\mathrm{\μ}}_s^{\mathrm{em}},{\mathrm{\μ}}_a^{\mathrm{em}})=T({\stackrel{\→}{r}}^{\′}\→\stackrel{\→}{r}|{\hat{s}}^{\′},{\mathrm{\μ}}_s^{\mathrm{em}},{\mathrm{\μ}}_a^{\mathrm{em}})C({\hat{s}}^{\′}\→\hat{s}|\stackrel{\→}{r},{\mathrm{\μ}}_s^{\mathrm{em}},{\mathrm{\μ}}_a^{\mathrm{em}})\\ =T({\stackrel{\→}{r}}^{\′}\→\stackrel{\→}{r}|{\hat{s}}^{\′},{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})\exp (-{\∫}_0^{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}({\mathrm{\μ}}_t^{\mathrm{em}}({\stackrel{\→}{r}}^{\′}+l\hat{s})-{\mathrm{\μ}}_t^{\mathrm{ex}}({\stackrel{\→}{r}}^{\′}+l\hat{s}))\mathrm{dl})\×\\ C({\hat{s}}^{\′}\→\hat{s}|\stackrel{\→}{r},{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})\frac{{\mathrm{\μ}}_s^{\mathrm{em}}\left(\stackrel{\→}{r}\right)}{{\mathrm{\μ}}_s^{\mathrm{ex}}\left(\stackrel{\→}{r}\right)}\\ =\frac{{\mathrm{\μ}}_s^{\mathrm{em}}\left(\stackrel{\→}{r}\right)}{{\mathrm{\μ}}_s^{\mathrm{ex}}\left(\stackrel{\→}{r}\right)}\exp (-{\∫}_0^{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}({\mathrm{\μ}}_t^{\mathrm{em}}({\stackrel{\→}{r}}^{\′}+l\hat{s})-{\mathrm{\μ}}_t^{\mathrm{ex}}({\stackrel{\→}{r}}^{\′}+l\hat{s}))\mathrm{dl})K(p^{\′}\→p,{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})\end{array}$

(4) throw in exciting light photon, follow the trail of photon transmission in biological tissues, utilize scattering coefficient to sample along exciting light path, the fluorescence intensity on calculating detector;

Following the trail of photon transmitting procedure concrete steps is in biological tissues:

(4.1) throw in exciting light photon, exciting light photon initial weight is set to 1.

(4.2) utilize the scattering coefficient of exciting light to sample all the time, the step-length that photon often walks is , scattering direction is determined by Henyey-Greenstein function.

(4.3) photon bag is decayed continuously along the scattering path of exciting light.

(4.4) exciting light photon is in non-fluorescence region, and the attenuation ratio of photon bag weight is exciting light photon is at fluorescence area, and the attenuation ratio of photon bag weight is

(4.5) according to relation between the kernel function of fluorescence excitation process and the kernel function of exciting light communication process, the probability of exciting light photon equilibrium state with the probability proportion of the fluorescent photon scattering excited is after fluorescent photon excites, still move along former exciting light scattering direction.

(4.6) if the optical parametric of fluorescence (absorption coefficient and scattering coefficient) is inconsistent with the optical parametric of exciting light, according to relation between the kernel function of fluorescence excitation process and the kernel function of exciting light communication process, the ratio that the decay of fluorescent photon photon bag weight in non-fluorescence region and exciting light attached bag weight decay is the ratio that fluorescent photon is decayed in decay and the exciting light attached bag weight of fluorescence area photon bag weight is $\exp \left((-{\mathrm{\μ}}_a^{\mathrm{ex}}+{\mathrm{\μ}}_{\mathrm{af}})l\right)\exp \left(({\mathrm{\μ}}_t^{\mathrm{em}}-{\mathrm{\μ}}_t^{\mathrm{ex}})l\right)\frac{{\mathrm{\μ}}_s^{\mathrm{em}}}{{\mathrm{\μ}}_s^{\mathrm{ex}}}.$

(4.7) along exciting light path, the fluorescence intensity on calculating detector.Through m collision rift, be y at the emission density function of p point _m(p):

$\begin{array}{c}{y}_m\left(p\right)=\underset{m=0}{\overset{\∞}{\mathrm{\Σ}}}\underset{i=0}{\overset{m}{\mathrm{\Σ}}}\∫...\∫S\left(p_0\right)K(p_0\→p_1,{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}}+{\mathrm{\μ}}_{\mathrm{af}})...\\ \×K(p_{i-1}\→p_i,{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}}+{\mathrm{\μ}}_{\mathrm{af}})\×K^{\mathrm{xm}}(p_i\→p_{i+1},{\mathrm{\μ}}_{\mathrm{af}})\×K(p_{i+1}\→p_{i+2},{\mathrm{\μ}}_s^{\mathrm{em}},{\mathrm{\μ}}_a^{\mathrm{em}})...\\ K(p_{m-1}\→p_m,{\mathrm{\μ}}_s^{\mathrm{em}},{\mathrm{\μ}}_a^{\mathrm{em}})\underset{k=0}{\overset{m}{\mathrm{\Π}}}{\mathrm{dp}}_k\end{array}$

Meanwhile, the probability density τ of photon can be obtained _mp () is expressed as:

${\mathrm{\τ}}_m\left(p\right)=S\left(p_0\right)K(p_0\→p_1,{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})...K(p_{m-1}\→p_m,{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})$

The weighting function w of photon _mp () is expressed as:

$\begin{array}{c}{w}_m\left(p\right)=g\left(p\right)\underset{i=0}{\overset{m}{\mathrm{\Σ}}}\exp (-{\∫}_0^{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}{\mathrm{\μ}}_{\mathrm{af}}({\stackrel{\→}{r}}^{\′}+l\hat{s})\mathrm{dl})\frac{{\mathrm{\η\μ}}_{\mathrm{af}}\left({\stackrel{\→}{r}}_i\right)}{{\mathrm{\μ}}_s^{\mathrm{ex}}\left({\stackrel{\→}{r}}_i\right)}\\ \×\frac{P_I({\hat{s}}_i\·{\hat{s}}_{i+1})}{P_A({\hat{s}}_i\·{\hat{s}}_{i+1})}\underset{j=i+1}{\overset{m}{\mathrm{\Π}}}\left(\frac{{\mathrm{\μ}}_s^{\mathrm{em}}\left({\stackrel{\→}{r}}_j\right)}{{\mathrm{\μ}}_s^{\mathrm{ex}}\left({\stackrel{\→}{r}}_j\right)}\right)\exp (-{\∫}_0^{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}({\mathrm{\μ}}_t^{\mathrm{em}}({\stackrel{\→}{r}}^{\′}+l\hat{s})-{\mathrm{\μ}}_t^{\mathrm{ex}}({\stackrel{\→}{r}}^{\′}+l\hat{s}))\mathrm{dl})\end{array}$

The function of fluorescence statistic of g (p) for detector detects.Fluorescence intensity on detector can calculate:

$\begin{array}{c}D={\∫y}_m\left(p\right)g\left(p\right)\mathrm{dp}\\ =\underset{m=0}{\overset{\∞}{\mathrm{\Σ}}}\∫...\∫{\mathrm{\τ}}_m\left(p\right)w_m\left(p\right)\underset{k=0}{\overset{m}{\mathrm{\Π}}}{\mathrm{dp}}_k\mathrm{dp}\end{array}$

The present invention is set forth further below by example.

Embodiment:

As shown in Figure 2, this cylinder model is full of the fat emulsion solution of 1% to cylinder model constructed in simulation, and the glass bar being full of fluorescent dye is inserted in fat emulsion solution.Carry out voxel subdivision to cylinder model, altogether containing 301401 voxels in the voxel model obtained after subdivision, each voxel is of a size of 0.05cm, and model is of a size of 3.05cm × 3.05cm × 4.05cm.The position of light source is taken at shown position, and detector is taken in the region on 180 degree, opposite, source, and it is evenly distributed on 80 layers, every layer of 180 detector.We list the tissue optical parameter value under exciting light and fluorescence bands in Table 1, and all the other parameters such as g is set to about 0.9, and refractive index is set to 1.37, the representative value of biological group under being near infrared spectrum.The photon number of simulation is 10 ⁹.

Tissue optical parameter value under table 1 exciting light and fluorescence bands

Fig. 3 a is normalization fluorescence intensity profile on the detector that obtains in the simulation of standard fluorescence Monte Carlo (sfMC) method, Fig. 3 b to separate in the simulation of coupling fluorescence Monte Carlo (dfMC) method normalization fluorescence intensity profile on the detector that obtains, and Fig. 4 is that the level line that standard conciliates normalization fluorescence intensity profile picture on the detector that obtains in the monte carlo modelling of coupling fluorescence compares.From detector, normalization fluorescence intensity profile and contour map can find out that two kinds of methods and resultses meet very well.

Claims (6)

1. separate a method for building up for coupling fluorescence monte-Carlo model, it is characterized in that comprising the following steps:

Step 1: build Particle transport equations and describe photon transmitting procedure in biological tissues under stable state;

Step 2: the emission density integral equation building transport equation equivalence describes process that exciting light photon transmits in biological tissues and exciting light photon is absorbed and inspires fluorescent photon, and fluorescent photon continues the process transmitted in biological tissues;

Step 3: the kernel function in fluorescence excitation process in the neumann series solution of integral equation and fluorescence communication process is all represented by the kernel function of exciting light communication process, thus solution is coupled out fluorescence excitation and the optical parametric of communication process and associating of path probability density function;

Step 4: throw in exciting light photon, follows the trail of photon transmission in biological tissues, utilizes exciting light scattering coefficient to sample, along exciting light path, and the fluorescence intensity on calculating detector.

2. the method for building up separating coupling fluorescence monte-Carlo model according to claim 1, is characterized in that in step 1, under stable state, photon transmitting procedure is in biological tissues described as:

$\begin{array}{c}\hat{s}\·\▿\mathrm{\φ}(\stackrel{\→}{r},\hat{s},v,t)+{\mathrm{\μ}}_t(\stackrel{\→}{r},v,t)\mathrm{\φ}(\stackrel{\→}{r},\hat{s},v,t)=S(\stackrel{\→}{r},\hat{s},v,t)+\\ {\∫}_{4\mathrm{\π}}{\mathrm{\μ}}_t(\stackrel{\→}{r},v,t)\mathrm{\φ}(\stackrel{\→}{r},{\hat{s}}^{\′},v,t)\·C({\hat{s}}^{\′}\→\hat{s})d{\mathrm{\Ω}}^{\′}\end{array}$

for the direction of photon current state;

for the direction of the next state of photon current state;

R is the position of photon current state;

V is the frequency of photon current state;

The solid angle unit of the next state that d Ω ' is photon current state;

for in t, near r point, in cell frequency interval and unit solid angle unit, direction exists neighbouring photon number;

for being photon in t, be emitted in the neighbouring cell frequency interval of r point and unit solid angle unit by external source, direction exists neighbouring photon number;

for direction is photon, at r point, place enters collision, and collision rift direction exists the photon probability number in direction;

μ _tfor extinction coefficient.

3. the method for building up separating coupling fluorescence monte-Carlo model according to claim 1, is characterized in that the process that in step 2, exciting light transmission and fluorescence excitation also transmit is described as respectively:

$x\left(p\right)=\∫x\left(p^{\′}\right)K(p^{\′}\→p,{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}}+{\mathrm{\μ}}_{\mathrm{af}}){\mathrm{dp}}^{\′}+S\left(p\right)$

$y\left(p\right)=\∫y\left(p^{\′}\right)K(p^{\′}\→p,{\mathrm{\μ}}_s^{\mathrm{em}},{\mathrm{\μ}}_a^{\mathrm{em}}){\mathrm{dp}}^{\′}+\∫x\left(p^{\′}\right)K(p^{\′}\→{p,\mathrm{\μ}}_{\mathrm{af}}){\mathrm{dp}}^{\′}$

P is photon transmission state, and p ' is the photon transmission state of after p;

The emission density that the emission density that x (p) is exciting light, y (p) are fluorescence;

S (p) is source emission density;

K is the kernel function of light communication process;

for exciting light scattering coefficient, for the absorption coefficient of exciting light;

for fluorescence scattering coefficient, for fluorescent absorption coefficient;

μ _affor the absorption coefficient of fluorophore.

4. the method for building up separating coupling fluorescence monte-Carlo model according to claim 1, to it is characterized in that in step 3 that solution is coupled out to excite and the fluorescent optics parameter of communication process and associating of path probability density function, the kernel function in fluorescence excitation process and fluorescence communication process all represents by the kernel function of exciting light communication process:

(1) kernel function of the kernel function exciting light communication process of fluorescence excitation process is expressed as:

$\begin{array}{c}K(p^{\′}\→p,{\mathrm{\μ}}_{\mathrm{af}})=T({\stackrel{\→}{r}}^{\′}\→\stackrel{\→}{r}|{\hat{s}}^{\′},{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}}+{\mathrm{\μ}}_{\mathrm{af}})C({\hat{s}}^{\′}\→\hat{s}|\stackrel{\→}{r},{\mathrm{\μ}}_{\mathrm{af}})\\ =K({\stackrel{\→}{r}}^{\′}\→\stackrel{\→}{r}|{\hat{s}}^{\′},{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})\exp (-{\∫}_0^{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}{\mathrm{\μ}}_{\mathrm{af}}({\stackrel{\→}{r}}^{\′}+l\hat{s})\mathrm{dl})\×\\ \frac{{\mathrm{\μ}}_t^{\mathrm{ex}}\left(\stackrel{\→}{r}\right)+{\mathrm{\μ}}_{\mathrm{af}}\left(\stackrel{\→}{r}\right)}{{\mathrm{\μ}}_r^{\mathrm{ex}}\left(\stackrel{\→}{r}\right)}\frac{C({\hat{s}}^{\′}\→\hat{s}|\stackrel{\→}{r},{\mathrm{\μ}}_{\mathrm{af}})}{C({\hat{s}}^{\′}\→\hat{s}|\stackrel{\→}{r},{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})}\\ =K(p^{\′}\→{p,\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})\exp (-{\∫}_0^{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}{\mathrm{\μ}}_{\mathrm{af}}({\stackrel{\→}{r}}^{\′}+l\hat{s})\mathrm{dl})\×\\ \frac{{\mathrm{\η\μ}}_{\mathrm{af}}\left(\stackrel{\→}{r}\right)}{{\mathrm{\μ}}_s^{\mathrm{ex}}\left(\stackrel{\→}{r}\right)}\frac{P_I({\hat{s}}^{\′}\·\hat{s})}{P_A({\hat{s}}^{\′}\·\hat{s})}\end{array}$

The position of the next state that r ' is photon current state;

---be the photon sent from r ', enter the density of collision at r place;

for the extinction coefficient of exciting light;

L is the step-length that photon often walks;

η is quantum efficiency;

P _ifor isotropic phase bit function, P _afor anisotropy phase function;

(2) kernel function of the kernel function exciting light communication process of fluorescence communication process is expressed as:

$\begin{array}{c}K(p^{\′}\→p,{\mathrm{\μ}}_s^{\mathrm{em}},{\mathrm{\μ}}_a^{\mathrm{em}})=T({\stackrel{\→}{r}}^{\′}\→\stackrel{\→}{r}|{\hat{s}}^{\′},{\mathrm{\μ}}_s^{\mathrm{em}},{\mathrm{\μ}}_a^{\mathrm{em}})C({\hat{s}}^{\′}\→\hat{s}|\stackrel{\→}{r},{\mathrm{\μ}}_s^{\mathrm{em}},{\mathrm{\μ}}_a^{\mathrm{em}})\\ =T({\stackrel{\→}{r}}^{\′}\→\stackrel{\→}{r}|{\hat{s}}^{\′},{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})\exp (-{\∫}_0^{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}({\mathrm{\μ}}_t^{\mathrm{em}}({\stackrel{\→}{r}}^{\′}+l\hat{s})-{\mathrm{\μ}}_t^{\mathrm{ex}}({\stackrel{\→}{r}}^{\′}+l\hat{s}))\mathrm{dl})\×\\ C({\hat{s}}^{\′}\→\hat{s}|\stackrel{\→}{r},{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})\frac{{\mathrm{\μ}}_s^{\mathrm{em}}\left(\stackrel{\→}{r}\right)}{{\mathrm{\μ}}_s^{\mathrm{ex}}\left(\stackrel{\→}{r}\right)}\\ =\frac{{\mathrm{\μ}}_s^{\mathrm{em}}\left(\stackrel{\→}{r}\right)}{{\mathrm{\μ}}_s^{\mathrm{ex}}\left(\stackrel{\→}{r}\right)}\exp (-{\∫}_0^{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}({\mathrm{\μ}}_t^{\mathrm{em}}({\stackrel{\→}{r}}^{\′}+l\hat{s})-{\mathrm{\μ}}_t^{\mathrm{ex}}({\stackrel{\→}{r}}^{\′}+l\hat{s}))\mathrm{dl})K(p^{\′}\→{p,\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})\end{array}$

for the extinction coefficient of fluorescence.

5. the method for building up separating coupling fluorescence monte-Carlo model according to claim 1, it is characterized in that the fluorescent photon path of sampling in step 4 is consistent with excitation light subpath, the weighting function of fluorescent photon is associated with the fluorescent optics parameter excited with communication process; Along excitation light subpath, directly calculate the weight of fluorescent photon according to fluorescence weighting function.

6. solution coupling fluorescence monte-Carlo model according to claim 1, is characterized in that the fluorescence intensity in step 4 on detector is expressed as:

$\begin{array}{c}D=\∫y_m\left(p\right)g\left(p\right)\mathrm{dp}\\ =\underset{m=0}{\overset{\∞}{\mathrm{\Σ}}}\∫...{\∫\mathrm{\τ}}_m\left(p\right)w_m\left(p\right)\underset{k=0}{\overset{m}{\mathrm{\Π}}}{\mathrm{dp}}_k\mathrm{dp}\end{array}$

D is the fluorescent photon intensity that detection receives;

The number of times that m collides for photon before arrival detector;

The function of fluorescence statistic of g (p) for detector detects;

τ _mp probability density that () is photon;

W _mp weighting function that () is photon;

(1) the probability density τ of photon _mp () is expressed as:

${\mathrm{\τ}}_m\left(p\right)=S\left(p_0\right)K(p_0\→p_1,{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})...K(p_{m-1}\→p_m,{\mathrm{\μ}}_s^{\mathrm{ex}},{\mathrm{\μ}}_a^{\mathrm{ex}})$

Wherein S (p)=S (p ₀), S (p ₀) be source emission density;

(2) the weighting function w of photon _mp () is expressed as:

$\begin{array}{c}{w}_m\left(p\right)=g\left(p\right)\underset{i=0}{\overset{m}{\mathrm{\Σ}}}\exp (-{\∫}_0^{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}{\mathrm{\μ}}_{\mathrm{af}}({\stackrel{\→}{r}}^{\′}+l\hat{s})\mathrm{dl})\frac{{\mathrm{\η\μ}}_{\mathrm{af}}\left(\stackrel{\→}{r_i}\right)}{{\mathrm{\μ}}_s^{\mathrm{ex}}\left(\stackrel{\→}{r_i}\right)}\\ \×\frac{P_I({\hat{s}}_i\·{\hat{s}}_{s+1})}{P_A({\hat{s}}_i\·{\hat{s}}_{i+1})}\underset{j=i+1}{\overset{m}{\mathrm{\Π}}}\left(\frac{{\mathrm{\μ}}_s^{\mathrm{em}}\left({\stackrel{\→}{r}}_j\right)}{{\mathrm{\μ}}_s^{\mathrm{ex}}\left({\stackrel{\→}{r}}_j\right)}\right)\exp \left({-\∫}_0^{|\stackrel{\→}{r}-{\stackrel{\→}{r}}^{\′}|}({\mathrm{\μ}}_t^{\mathrm{em}}({\stackrel{\→}{r}}^{\′}+l\hat{s})-{\mathrm{\μ}}_t^{\mathrm{ex}}({\stackrel{\→}{r}}^{\′}+l\hat{s}))\mathrm{dl}\right)\end{array}$