Message Boards Message Boards


Problem: Light transport through biological tissue

Posted 1 year ago
6 Replies
0 Total Likes

I'm new to Mathematica and while it seems like there should be a way to do this (it is effectively simultaneous equations) I don't know how to input it to Mathematica. I've tried using solve and setting it in a table first but it seems I'm not realizing something or missing something in the documentation about how to play with a few.

The problem represents light passing through tissue. There are two variables in the equation that I need to know, a and b. (if there is a method to input equations/code here can someone point it out to me?)

enter image description here

I(lam), I0(lam) are measured values. I will know these, they are the intensity of the light coming out and going into the volume at wavelength lambda. M(lam) is the absorption of tissue at the wavelength lamda, this I will also know at all lambdas.

d thickness of the tissue, this is a set value.

a and b are what I would like to know. That term represents the light scattering function for the tissue.

What I would like mathematica to do is tell me how many lambda's, or colors of light, do I have to measure at to find what a and b are, and the equation that gives me a and b.

6 Replies

If it's measured data, you probably want one of the *Fit family of function, for example FindFit or NonlinearModelFit. How many values are needed should depend on how noisy the data is, but statistics is not a field I am familiar enough with to give a good answer. The FittedModel returned by NonlinearModelFit will report various measurements of quality; see the docs.

Hi Alan

Yours problem can be solved if you have 2 equations.

Let's try:

Solve[a + b == c, {a, b}]

(* Solve::svars: Equations may not give solutions for all "solve" variables.*)

Give an error occurs because the number of variables is greater than the number of equations.

If I have 2 equations then:

 Solve[{a + b == c, a == 1}, {a, b}]

 (* {{a -> 1, b -> -1 + c}}*)

is Ok.


Welcome to Wolfram Community! Please make sure you know the rules:

Please do NOT use images for posting code.

The rules explain how to format your code properly. If you do not format code, it may become corrupted and useless to other members. Please EDIT your posts and make sure code blocks start on a new paragraph and look framed and colored like this.

int = Integrate[1/(x^3 - 1), x];
Map[Framed, int, Infinity]

enter image description here

You can do it a bit more directly. The idea is to get an expression in a and b which could be evaluated without major problems.

Start with your equation ( I write w (wavelength) for lambda to avoid the \ [...] characters)

eq1 = i[w] == i0[w] Exp[ - m[w] d/2 ((3 a w^-b)/m[w])^(1/2)]

Now I apply certain transforms of Logarithms to each side of eq1 (that is done by f/@eq1. Try a simple example)

eq2 = (Log[#] /. Log[a_ b_] -> Log[a] + Log[b] /. Log[Exp[x_]] -> x) & /@ eq1

and another one

eq3 = - (# - Log[i0[w]]) & /@ eq2 /. Log[x] - Log[y] -> Log[x/y]

to get as eq3

Log[i0[w]/i[w]] == 1/2 Sqrt[3] d Sqrt[(a w^-b)/m[w]] m[w]

Now take again logarithms on both sides

eq4 = Log[#] & /@ eq3


Log[Log[i0[w]/i[w]]] == Log[1/2 Sqrt[3] d Sqrt[(a w^-b)/m[w]] m[w]]

To simplify or expand the expression on the right hand side I define a function

flog[Log[a_]] := Plus @@ (Log /@ List @@ a)

Look what it does

In[7]:= flog[Log[1/2 Sqrt[3] d Sqrt[(a w^-b)/m[w]] m[w]]]

Out[7]= -Log[2] + Log[3]/2 + Log[d] + Log[Sqrt[(a w^-b)/m[w]]] +  Log[m[w]]

Apply it to the right side of eq 4 (try what MapAt[h, eq4, 2] does) and simplify further

eq5 = MapAt[flog, eq4, 2] /. Log[Sqrt[x_]] -> 1/2 Log[x]

Isolate the terms with a and b on the right side

eq6 = 2 (# + Log[2] - Log[3]/2 - Log[m[w]] - Log[d]) & /@ eq5

eq7 = MapAt[flog, eq6, 2] /. Log[a_^x_] -> x Log[a]

eq8 = (# + Log[m[w]]) & /@ eq7

resulting in

-2 Log[d] + 2 Log[Log[i0[w]/i[w]]] - Log[(3 m[w])/4] == Log[a] - b Log[w]

(I am quite aware that all this could have been done "by hand"). Anyhow, in this expression the left hand side ( which could - of course - be simplified further) is known as it contains known values only, the right hand has the unknowns b and a in the form of Log[ a ]. So with a couple of measurements you should be able to find the values you are looking for, preferably by a least square fit.

To simplify the lefthand side of eq8 use

eq8[[1]] /. -2 Log[d] -> 2 Log[1/d] /. x_ Log[a_] -> Log[a^x] //. Log[a_] + Log[b_] -> Log[ a b]

With two unknowns a and b you might expect that two measurements (at two colors or wavelengths) will give your desired answer.

Simplifying the left hand side of the eq8 this will do for two measurements

t1 = Log[(4 Log[i0[w]/i[w]]^2)/(3 d^2 m[w])] /. w -> w1;
t2 = Log[(4 Log[i0[w]/i[w]]^2)/(3 d^2 m[w])] /. w -> w2;
Solve[{t1 == Log[a] - b Log[w1], t2 == Log[a] - b Log[w2]}, {a, b}] // FullSimplify

The resulting expression is somewhat complicated, so you should better feed numbers in.

But I am afraid just two measurement will not be sufficient do determine a and b satisfactorily. So I think a couple of measurements and a list square fit are necessary.

Posted 1 year ago

Thank you to everyone for your help with this!

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract