75048 Mathematical modeling 1998

Lin & Segel’s physiological flow problem

Harald Hanche-Olsen

                       hanche@math.ntnu.no

A brief description of the model

Osmosis is the transport of water through a membrane resulting from differences in the solute concentrations, C_{1} vs. C_{2}, of the two sides of the membrane. The resulting flow rate, in volume per unit area and time unit is given by

J = P (C_{2} − C_{1})  (1)

Here the units are: [J] = m^{3} · m^{−2}s^{−1} = m s^{−1}; [C_{i}] = osmol m^{−3}; and hence [P] = m^{4}osmol^{−1}s^{−1}. Diffusion is given by Fick’s law: The solute flux is proportional to the concentration gradient,

F_{D} = −D \frac{∂C}{∂x}  (2)

where units are [F_{D}] = osmol m^{−2}s^{−1}; and [D] = m^{2}s^{−1}. D is called the diffusion coefficient. Chemical pumps at the inner end of the channel move solute into the channel via a chemical pump whose mechanism we will not try to understand here. This is postulated to take place at a constant rate N_{0} per area and time, that is, [N_{0}] = osmol m^{−2}s^{−1}.

We model the channel as a long thin channel of constant crosssectional area a, circumference c, and length L (see figure 1). The chemical pump acts only at the inner end of the channel, in a section of length δ. Since the channel is assumed to be thin, we will assume all quantities to be constant over any cross section of the channel, so that we need only a single space coordinate, namely x^{∗} (with x^{∗} = 0 at the innermost end, and x^{∗} = L at the outlet).

We look for a stationary state for the system, so that all unknown quantities are functions of x^{∗} alone. Unknowns are the solute concentration C^{∗} and the water velocity v^{∗}.

[figure: δ; L; x^{∗} = 0; δ; x^{∗} = L]

Figure 1: The narrow channel

Another quantity of interest is the solute concentration C_{0} in the surrounding tissue, which is also assumed to be the solute concentration at the outlet of the channel.

Mass balance for the water. Since water is neither created nor destroyed, the net amount of water leaving the channel between x_{1}^{∗} and x_{2}^{∗} (figure 2) must be zero:

av^{∗}(x_{2}^{∗}) − av^{∗}(x_{1}^{∗}) − ∫_{x_{1}^{∗}}^{x_{2}^{∗}} cJ(x^{∗})dx^{∗} = 0

where J is given by (1) (with C_{2} = C^{∗} and C_{1} = C_{0}). This can be

[figure: x_{1}^{∗}; x_{2}^{∗}; av^{∗}(x_{1}^{∗}); av^{∗}(x_{2}^{∗}); ∫_{x_{1}^{∗}}^{x_{2}^{∗}} P (C^{∗} − C_{0})c dx^{∗}]

Figure 2: Conservation of water rewritten using the fundamental theorem of calculus into

∫_{x_{1}^{∗}}^{x_{2}^{∗}} (a\frac{dv^{∗}}{dx^{∗}} − cJ(x^{∗})) dx^{∗} = 0

and so the integrand must be zero if this is to hold for every choice of x_{1}^{∗} and x_{2}^{∗}:

\frac{dv^{∗}}{dx^{∗}} = \frac{c}{a}J(x^{∗}) = \frac{c}{a}P (C^{∗} − C_{0})  (3)

Mass balance for solute. Solute, like water, is neither created nor destroyed in this process. Figure 3 applies the the innermost region x^{∗} < δ; when x^{∗} > δ on the other hand, the contribution from the term N_{0}c(x_{2}^{∗} − x_{1}^{∗}) should be replaced by 0. F^{∗} represents the total

[figure: x_{1}^{∗}; x_{2}^{∗}; aF^{∗}(x_{1}^{∗}); aF^{∗}(x_{2}^{∗}); N_{0}c(x_{2}^{∗} − x_{1}^{∗})]

Figure 3: Conservation of solute

solute flux in the channel, and is the sum of a diffusive and a convective term: F^{∗} = F_{D}^{∗} +F_{K}^{∗} where F_{D}^{∗} is given by (2) and F_{K}^{∗} = v^{∗}C^{∗}. A quick glance at the figure leads to the equation

aF^{∗}(x_{2}^{∗}) − aF^{∗}(x_{1}^{∗}) − N_{0}c(x_{2}^{∗} − x_{1}^{∗}) = 0, 0 ≤ x_{1}^{∗} < x_{2}^{∗} < δ

We could of course change this to a differential equation as above, but instead we can use it to find an expression for F^{∗}(x^{∗}) by setting x_{1}^{∗} = 0 and x_{2}^{∗} = x^{∗}. Because F^{∗}(0) = 0 we get aF^{∗}(x^{∗}) = cN_{0}x^{∗}. For δ ≤ x_{1}^{∗} < x_{2}^{∗} < L it is simpler still, since the term N_{0}c(x_{2}^{∗} − x_{1}^{∗}) is dropped, and we conclude that F^{∗}(x^{∗}) is constant for x^{∗} ≥ δ. Thus we find

v^{∗}C^{∗} − D \frac{dC^{∗}}{dx^{∗}} = F^{∗}(x^{∗}) = \begin{cases}\frac{c}{a}N_{0}x^{∗} & for 0 ≤ x^{∗} ≤ δ \\ \frac{c}{a}N_{0}δ & for δ ≤ x^{∗} ≤ L \end{cases}  (4)

Boundary conditions. In deriving (4) we already used one boundary condition, namely F^{∗}(0) = 0. Since this is automatically fulfilled when (4) holds, we no longer need consider this. Thus we are left with two boundary conditions (5) and (6), plus the interior boundary conditions (7):

\begin{aligned}v^{∗}(0) & = 0 & (5) \\ C^{∗}(L) & = C_{0} & (6) \end{aligned}

v^{∗}, C^{∗}, \frac{dC^{∗}}{dx^{∗}} continuous at x^{∗} = δ  (7)

The inner boundary conditions need to be added because (4) introduces a potential discontinuity into the problem for x^{∗} = δ. The continuity of v^{∗} is necessary for water conservation, and F^{∗} = v^{∗}C^{∗} − D(dC^{∗}/dx^{∗}) must be continuous for the sake of solute conservation. C^{∗} must be continuous because diffusion will make it so, and it follows that dC^{∗}/dx^{∗} is continuous as well.

Dimensional analysis and scaling

We have arrived at a model given by equations (3), (4), (5), (6), (7). The parameters of the problem occur only in the combinations cP/a, cN_{0}/a, C_{0}, D, δ and L. which does simplify dimensional analysis somewhat. We start with these parameter combinations and their units:

δ  C_{0}  \frac{c}{a}N_{0}  \frac{c}{a}P  D L ; m \frac{osmol}{m^{3}} \frac{osmol}{m^{3}s} \frac{m^{3}}{osmol s} \frac{m^{2}}{s} m

From these six quantities we can, as expected, produce three independent dimensionless combinations, for example

ν = \frac{N_{0}}{P C_{0}^{2}} , η = \frac{aC_{0}D}{cN_{0}δ^{2}} , λ = \frac{L}{δ}

Scaling. The solute flux must clearly change from 0 to its maximum value in the region from x^{∗} = 0 to x^{∗} = δ. This ought to make δ a good length scale. We need a suitable scale for v^{∗}. Unless the diffusion term dominates completely in (4), the other two terms in that equation must be the same order of magnitude, which we write v^{∗}C^{∗} ∼ cN_{0}δ/a. Assuming C^{∗} ∼ C_{0} we then get v^{∗} ∼ cN_{0}δ/(aC_{0}), which we chose as a scale for v^{∗}. So we scale the problem using

C^{∗} = C_{0}C
v^{∗} = \frac{cN_{0}δ}{aC_{0}} v  (8)
x^{∗} = δx

and the complete problem becomes (from now on writing C′ = dC/dx, etc)

\begin{aligned}νv′ & = C − 1 \\ Cv − ηC′ & = \begin{cases}x & for 0 ≤ x ≤ 1 \\ 1 & for 1 ≤ x ≤ λ \end{cases} \\ v(0) & = 0 \\ C(λ) & = 1 \end{aligned}

v, C, C′ continuous at x = 1

Perturbation analysis

A somewhat naïve attempt at a perturbation analysis is to treat ν as a small parameter, and to express v and C as power series in ν. If we put v = v_{0} + νv_{1} + … and C = G_{0} + νG_{1} + …,^{1} substitute these in the model, and compute to second order in ν we find

νv′_{0} = G_{0} − 1 + νG_{1} + 𝒪(ν^{2}) ⟹ G_{0} = 1, v′_{0} = G_{1}

1The symbol C_{0} is already in use, so I use G_{i} for the coefficients of the power series for C. The alternative is Lin & Segel’s C^{(i)}, which I find overly complex.

and

G_{0}v_{0} + ν(G_{1}v_{0} + G_{0}v_{1}) − ηG′_{0} − νηG′_{1} = \begin{cases}x + 𝒪(ν^{2}) & for 0 ≤ x ≤ 1 \\ 1 + 𝒪(ν^{2}) & for 1 ≤ x ≤ λ \end{cases}

from which v_{0} = x for x < 1 and v_{0} = 1 for x > 1, and hence G_{1} = 1 for x < 1 and G_{1} = 0 for x > 1. But this contradicts the requirement of continuity of C at x = 1!

A more successful attempt. From a study of the numerical values of the parameters we see that N_{0} is the most variable of all the quantities in the problem, ranging over several orders of magnitude. If we change N_{0} so that ν becomes small we see that η has a tendency to become large at the same time, with νη approximately constant. It might, therefore, be more useful to assume νη constant in the perturbation expansion. As it turns out, the final answer becomes a bit simplify if we name this constant as follows:

νη = ψ^{−2}

Elimination of η from the problem now converts it into

\begin{aligned}ν \frac{dv}{dx} & = C − 1 \\ νCv − \frac{1}{ψ^{2}} \frac{dC}{dx} & = ν \begin{cases}x & for 0 ≤ x ≤ 1 \\ 1 & for 1 ≤ x ≤ λ \end{cases} \\ v(0) & = 0 \\ C(λ) & = 1 \end{aligned}

v, C, C′ continuous at x = 1

Again we put v = v_{0} + νv_{1} + ⋯ and C = G_{0} + νG_{1} + ⋯, and find

νv′_{0} = G_{0} − 1 + νG_{1} + 𝒪(ν^{2})
νG_{0}v_{0} − ψ^{−2}G′_{0} − νψ^{−2}G′_{1} = ν \begin{cases}x \\ 1 \end{cases}\} + 𝒪(ν^{2})
v_{0}(0) + νv_{1}(0) + ⋯ = 0
G_{0}(λ) + νG_{1}(λ) + ⋯ = 1

where the expression in braces is short for “pick the upper expression for x < 1, and the lower one for x > 1.” For the constant term, we find the two equations G_{0} − 1 = 0 and −ψ^{−2}G′_{0} = 0, which fortunately are consistent, with the solution G_{0} = 1. This also fits the boundary condition G_{0}(λ) = 1. For the first order term (coefficient of ν^{1}) we get

v′_{0} = G_{1}, v_{0}C_{0} − ψ^{−2}C′_{1} = \begin{cases}x \\ 1 \end{cases}\}

with the solution

v_{0} − ψ^{−2}v′′_{0} = \begin{cases}x \\ 1 \end{cases}\}

The general solution in 0 < x < 1 is v_{0} = x + Ae^{ψx} + Be^{−ψx}, but from v_{0}(0) = 0 we then get

v_{0}(x) = x − K_{1} sinh(ψx), 0 < x < 1

The general solution in 1 < x < λ is v_{0} = 1 + Ae^{ψx} + Be^{−ψx}, but since v′_{0} = G_{1} and G_{1}(λ) = 0 we find v′_{0}(λ) = 0, and we can write

v_{0}(x) = 1 − K_{2} cosh(ψ(x − λ)), 1 < x < λ

Since v′_{0} = G_{1} and C is continuous at x = 1, we see that not only is v_{0} continuous at x = 1, but so is v′_{0}. This produces two equations for the determination of K_{1} and K_{2}:

1 − K_{1} sinh(ψ) = 1 − K_{2} cosh(ψ(1 − λ))
1 − K_{1}ψ cosh(ψ) = −K_{2}ψ sinh(ψ(1 − λ))

Rewrite this in matrix form and use the identity cosh(α−β) = cosh(α) cosh(β)− sinh(α) sinh(β) to discover the matrix has determinant cosh(ψλ), from which we quickly find the solution

(\begin{matrix}K_{1} \\ K_{2} \end{matrix}) = \frac{1}{ψ cosh(ψλ)} (\begin{matrix}cosh(ψ(1 − λ)) \\ sinh(ψ) \end{matrix})

using Cramer’s rule. If you wish to compare this with the results of Lin & Segel, you should note that their κ is my ψλ.

We are ultimately interested in Os^{∗}, which is defined as the total solute flux at the outlet divided by the total water flux there. Going back to equation (4), we see that F^{∗}(L) = (c/a)N_{0}δ, so that

Os^{∗} = \frac{cN_{0}δ}{av^{∗}(L)} = \frac{C_{0}}{v(λ)}

where we used (8). This might be a good point at which to introduce dimensionless osmolarity

Os = \frac{Os^{∗}}{C_{0}} = \frac{1}{v(λ)} = \frac{1}{1 − K_{2}} + 𝒪(ν)

In many cases of interest, ψ is relatively small, so that sinh(ψ) ≈ ψ. If this is introduced into the expression for K_{2} we get K_{2} ≈ 1/ cosh(ψλ), so that

Os ≈ \frac{cosh(κ)}{cosh(κ) − 1}

where κ = ψλ. For a discussion of these results and a comparison with numerical and experimental results I refer to Lin & Segel.

Reference

C.C Lin & L.A. Segel: Mathematics applied to deterministic problems in the natural sciences. ISBN 0898712297. Chapter 8 (pp. 244– 276).