Difusión de Fick e introducción a difusión en

sistemas confinados.
Reporte 1.
Roberto Verdel Aranda
27 de mayo de 2014

1. Planteamiento del problema

En términos generales, llamamos movimiento browniano 1 al movimiento
seguido por una partı́cula muy pequeña inmersa en un fluido. Este movi-
miento se caracteriza por ser continuo e irregular, y es provocado por los
choques entre la partı́cula y las moléculas que constituyen al fluido. Tal es el
movimiento que siguen muchos de los componentes celulares, siendo agua el
principal constituyente del medio en este caso. Cuando con este movimiento
se tiende a uniformizar la concentración de las de las partı́culas brownianas,
hablaremos de difusión. En general, la solución al problema de difusión, con-
siste en dar la concentración de las partı́culas como función de su posición y
del tiempo.
En esta primera parte del proyecto, discutiremos cómo resolver el proble-
ma de difusión en sistemas simples usando la ecuación de Fick. Abordaremos
también el problema de difusión en sistemas confinados en geometrı́as com-
plejas, donde presentaremos el método de los propagadores. En ambos casos,
discutiremos algunas aplicaciones en sistemas biológicos, aunque el trata-
miento presentado en este trabajo es general y puede emplearse de igual
modo en el estudio de otros fenómenos.
1
Llamamado ası́ en reconocimiento al botánico escocés Robert Brown (1773-1858),
quien en 1828 presentó el primer estudio detallado del fenómeno.

1
2. Teorı́a microscópica
Empezaremos el estudio de la difusión haciendo una descripción mi-
croscópica de este fenómeno. A manera de de nota histórica, cabe mencionar
que fueron Marian Smoluchowski (1872-1917) y Albert Einstein (1879-1955),
quienes a inicios del siglo pasado publicaron algunos trabajos dando una des-
cripción cualitativa y, por primera vez, cuantitativa del fenómeno en cuestión,
del tipo de la que daremos en seguida. Ambos cientı́ficos se basaron en ideas
de la teorı́a cinética de Maxwell y Boltzmann y llegaron esencialmente a
los mismos resultados, aunque por la manera en que fueron presentados, los
trabajos de Einstein tuvieron una mayor influencia que los de Smoluchowski.
Con vistas a favorecer la claridad en nuestra exposición, comenzaremos
considerando el problema en una dimensión, generalizando más adelante a
tres dimensiones. Por principio de cuentas, modelaremos la posición de una
partı́cula browniana usando caminatas aleatorias simples. El tiempo se con-
sidera discreto. Las hipótesis concomitantes a este modelo son las que se
enumeran a continuación:
1. Cada partı́cula browniana puede desplazarse a la izquierda o derecha
una distancia ∆x en cada paso de tiempo de duración ∆t. Supondremos
que ∆x y ∆t se mantienen fijas, por simplicidad.
2. Las partı́culas avanzan a la izquierda o derecha con la misma probabi-
lidad p = 1/2.
3. Aceptamos como válida la propiedad de Markov, i.e. que las partı́culas
no guardan memoria de lo que hicieron en el paso de tiempo anterior.
4. El movimiento de las partı́culas brownianas es debido a las colisiones
con las moléculas del fluido. Las partı́culas brownianas no interactúan
entre sı́, o sea que se mueven de manera independiente unas de otras.
(Hipótesis válida sólo en el caso de bajas densidades.)
Bajo estos supuestos, considérese un conjunto de N partı́culas brownia-
nas. Denotamos por xi (n) la posición de la i-ésima partı́cula después de n
pasos temporales. De acuerdo con la primera hipótesis, se tiene que,
xi (n) = xi (n − 1) ± ∆x, (2.1)
esto es, que la posición al paso n difiere de la posición en el paso anterior por
una cantidad ±∆x. El signo ± es debido a que la partı́cula puede desplazarse
a la derecha o izquierda.

2
 De las otras tres hipótesis es claro que después de varios pasos temporales,
en promedio, la mitad de las partı́culas se habrán movido una longitud +∆x,
y la otra mitad una distancia −∆x.
Tomemos el promedio2 de los desplazamientos de las N partı́culas después
de n pasos temporales:
N
1 X
hx(n)i = xi (n). (2.2)
N i=1

Introduciendo (2.1) en la ec. (2.2), y por lo dicho lı́neas arriba, tenemos

que,
N
1 X
hx(n)i = xi (n − 1) = hx(n − 1)i. (2.3)
N i=1
La ec. (2.3) nos dice que el desplazamiento promedio permanece inalte-
rado en los sucesivos pasos de tiempo. La única información que se extrae de
lo anterior, es que la distribución espacial de las partı́culas es simétrica con
respecto a su posición inicial para todo paso de tiempo. En particular, si al
inicio todas las partı́culas se encuentran en el origen, x = 0, después de n
pasos, se tiene que, hx(n)i = 0.
Necesitamos pues, otra cantidad que nos permita cuantificar de cierta
manera cuánto se esparcen las partı́culas al transcurrir el tiempo, respecto de
donde iniciaron. La cantidad buscada es el desplazamiento cudrático medio,
hx2 (n)i. Elevando al cuadrádo la ecuación (2.1) y sustituyendo el resultado
de ello en la definición (2.2), uno obtiene lo siguiente,
N
2 1 X 2
hx (n)i = x (n)
N i=1 i
N
1 X 2
= [xi (n − 1) + (∆x)2 ± 2∆x · xi (n − 1)]. (2.4)
N i=1

Nuevamente, como el promedio, h±∆xi = 0 (en promedio, la mitad de

las partı́culas se mueven a la derecha y la otra mitad a la izquierda), y con la
definición (2.2), la expresión anterior puede escribirse de forma más compacta
como,
hx2 (n)i = hx2 (n − 1)i + (∆x)2 . (2.5)
2
En todo el escrito, los paréntesis angulares denotarán un promedio.

3
 Poniendo n = 1 en la ecuación anterior, y suponiendo que las partı́cu-
las se encuentran inicialmente en el origen, se halla que el desplazamiento
cuadrático medio es (∆x)2 . Después de iterar n pasos, llegamos a que,

hx2 (n)i = n(∆x)2 . (2.6)

Como asignamos una duración ∆t a cada paso temporal, el tiempo total

transcurrido después de n pasos es, t = n∆t, de donde n = t/∆t, y la ecuación
(2.6) queda de la siguiente manera,

2 (∆x)2
hx (t)i = t. (2.7)
∆t
Definimos ahora el coeficiente de difusión:

(∆x)2
D := lı́m , (2.8)
∆x,∆t→0 2∆t

con lo que llegamos al siguiente resultado importante:

hx2 (t)i = 2Dt. (2.9)

Esta ecuación nos indica que el desplazamiento cuadrático medio se in-

crementa con el tiempo de forma proporcional a 2D. Una consecuencia im-
portante de este resultado es que a la partı́cula le toma poco tiempo recorrer
distancias cortas, pero le lleva una gran cantidad de tiempo hacer un reco-
rrido largo.

2.1. Ecuación de difusión de Fick

Nuestro siguiente paso es hacer la transición al continuo en el modelo
anterior. La estrategia a seguir es la siguiente: escribiremos una ecuación de
evolución para la función de probabilidad, Pn (j), de encontrar a la partı́cula
browniana en determinada posición en un momento dado. Luego, haremos
tender a cero las cantidades, ∆x y ∆t, para posteriormente realizar un desa-
rrollo en serie de Taylor en torno al cero.
Seguimos pues, con la partácula browniana restringida a moverse en una
dimensión, pudiendo desplazarse a uno de sus puntos vecinos más cercanos
en cada paso de tiempo. Denotamos por p la probabilidad de que la partı́cula
efectúe un movimiento a la derecha, y por q la probabilidad de que lo haga

4
hacia la izquierda. Estas probabilidades no necesariamente tienen que ser
iguales, pero deben de satisfacer que, p + q = 1 (la partı́cula se desplaza en
cada paso de tiempo). Si la partı́cula se encuentra en el paso n en el punto
etiquetado con el marcador j−1, entonces en el paso siguiente puede terminar
en la posición j-ésima con probabilidad p. Análogamente, si se encontraba
en j + 1, puede alcanzar el j-ésimo punto con probabilidad q. En definitiva,
tenemos que la probabilidad de hallar a la partı́cula en el punto j-ésimo al
paso n + 1, está dada por,

Pn+1 (j) = pPn (j − 1) + qPn (j + 1). (2.10)

Ahora pasaremos de los enteros al lı́mite continuo. Esto se hace sustitu-

yendo las variables j y n, por sus contrapartes continuas, x y t. Considerando
las relaciones: x = j∆x y t = n∆t, la ecuación (2.10) puede reescribirse como,

P (j∆x, (n + 1)∆t) = pP ((j − 1)∆x, n∆t) + qP ((j + 1)∆x, n∆t),

o bien,
P (x, t + ∆t) = pP (x − ∆x, t) + qP (x + ∆x, t). (2.11)
El lı́mite continuo se obtiene al hacer que tanto ∆x como ∆t, se aproximen
a cero. Considerando dicho lı́mite, se hace un desarrollo en serie de Taylor en
torno al cero de las funciones que aparecen en (2.11), con lo que se obtiene,

∂P (x, t)
P (x, t + ∆t) ≈ P (x, t) + ∆t + ··· ,
∂t
y,
∂P (x, t) (∆x)2 ∂ 2 P (x, t)
P (x ± ∆x, t) ≈ P (x, t) ± ∆x + + ··· .
∂x 2 ∂x2
Reemplazando estas expresiones en (2.11), se tiene que,

∂P (x, t) ∂P (x, t) (∆x)2 ∂ 2 P (x, t)

P (x, t) + ∆t = P (x, t) + (q − p)∆x + +··· ,
∂t ∂x 2 ∂x2
ya que p + q = 1. Desestimando términos de orden superior en los desarrollos
anteriores, y con la definiciones (2.8) y (2.12),

∆x
ν := − lı́m (q − p) , (2.12)
∆x,∆t→0 ∆t

5
uno llega a la expresión conocida como ecución de difusión con arrastre:

∂P (x, t) ∂ 2 P (x, t) ∂P (x, t)

=D 2
−ν . (2.13)
∂t ∂x ∂x
El arrastre es debido al segundo término en el miembro derecho de (2.13),
y es consecuencia de tener una dirección preferencial. Éste desaparece cuando
p = q = 1/2, ya que entonces ν = 0. En este último caso, la ecuación (2.13)
se reduce a la ecución de difusión de Fick 3 :

∂P (x, t) ∂ 2 P (x, t)
=D . (2.14)
∂t ∂x2

Consideremos ahora el caso tridimensional pensando que el fluido en que

se mueven las partı́culas brownianas es isótropo (con lo que las constantes
de difusión en las tres direcciones son iguales). Es inmediato ver que (2.14)
se transforma en
∂P (~r, t)
= D∇2 P (~r, t), (2.15)
∂t
donde ~r es el vector con componentes, x, y, z, y ∇2 es el operador laplaciano,
dado por
∂2 ∂2 ∂2
∇2 := + + . (2.16)
∂x2 ∂y 2 ∂z 2
En la constante de difusión que aparece en (2.15) ha de cambiarse el 2
del denominador de la definición (2.8) por un 6, ya ahora aparece sumado
tres veces el factor 2Dt en la relación para desplazamiento cuadrático medio
(siempre y cuando se tenga la misma constante de difusión en cada una de
las direcciones).
La solución a la ecuación (2.15) da la probabilidad de encontrar a una
partı́cula browniana en la posición especı́ficada por el vector ~r, al tiempo t.
Si consideramos nuevamente un conjunto de N partı́culas brownianas y mul-
tiplicamos la función P (~r, t) por el número N , obtendremos la concentración
de las partı́culas como una función del punto y del tiempo,

C(~r, t) := N P (~r, t). (2.17)

3
Derivada por primera vez en 1855 por el médico y fisiólogo alemán Adolf E. Fick
(1829-1901).

6
 Esta función describe la distribución de partı́culas en cualquier posición
a todo tiempo. Claramente, también satisface la ecuación de Fick.

∂C(~r, t)
= D∇2 C(~r, t). (2.18)
∂t

3. Teorı́a macroscópica
Daremos ahora una descricpción fenomenológica de la difusión. Comenza-
remos repasando conceptos tales como la conservación de partı́culas (masa),
flujo y corriente. Después deduciremos la primera ley de Fick, ası́ como la
ecuación de continuidad.
El primer concepto que debemos tener presente es el último que se intro-
dujo en la sección anterior, a saber, el concepto de concentración. La concen-
tración C(~r, t), pude definirse también como aquella función cuya integral de
volumen sobre un volumen V es igual al número de partı́culas contenidas en
dicho volumen, esto es,
Z
N (t) = d3~r C(~r, t). (3.1)
V

Se presenta ahora el concepto de flujo de partı́culas por unidad de área y

por unidad de tiempo, el cual está dado por

J~ = ~v C(V, t), (3.2)

donde ~v es la velocidad del flujo.

Para entender la expresión anterior pensemos en lo siguiente: un número
de partı́culas N que ocupan un volumen V , atraviesan un área A en un tiempo
t cuando se tiene que V = −A~n·~v t, con ~n el vector unitario perpendicular a A,
que apunta hacia el exterior. El signo menos se debe al hecho de que cuando
~v y ~n tienen componentes apuntando en direcciones opuestas, el número de
partı́culas crece en el seno del volumen. Si por ejemplo, ~v y ~n, son vectores
colineales, entonces ~n · ~v = ±v, donde se toma el signo + si los vectores
apunta en el mismo sentido, o el signo − en caso contrario. Sustituyendo este
volumen en (3.2), obtenemos,

~v N
−(±) ,
v At

7
donde el vector ~v /v solo da la dirección del flujo. De manera que se obtiene el
número de particulas que pasa por el área A en un tiempo t. Si las partı́culas
se mueven en la dirección en que apunta ~n, la expresión anterior lleva un signo
menos, lo cual, por la elección de ~n, significa que las parı́culas abandonan el
volumen considerado y viceversa.
Deduciremos ahora la primera ley de Fick, que muestra otra forma por
medio de cual están relacionados el flujo y la concentración. Para esto, con-
sidérese un recipiente cilı́ndrico de longitud 2L. Supondremos que el eje del
cilindro coincide con el eje x y que los extremos del mismo se encuentran en
x = ±L. Supongamos que divimos tal volumen en dos compartimentos igua-
les de longitud L. Estos volúmenes, digamos V1 y V2 , contienen un número
de partı́culas N1 y N2 , respectivamente. Estudiemos el cambio en el número
de partı́culas por unidad de área en el compartimento 1. Teniendo en cuenta
la ec. (3.1), se tiene que,
Z
1 dN1 1 d
= d3~r C(~r, t).
A dt A dt V1

Aquı́, A es el área de la circunferencia que delimita a los volúmenes. Supon-

dremos además que se tiene simetrı́a radial, con lo cual, la concentración sólo
depende de la coordenada x. En tal caso, la expresión anterior queda como
sigue,
1 d 0
Z Z 0
1 dN1 ∂C(x, t)
= dx A C(x, t) = dx ,
A dt A dt −L −L ∂t
donde en el último paso hemos podido introducir la derivada en el integrando
ya que los lı́mites de integración no dependen de t. Usando ahora la ecuación
de difusión [ec. (2.14)], en la expresión anterior, obtenemos que,
Z 0
∂ 2 C(x, t)
J= dx D
−L ∂x2
    
∂C(x, t) ∂C(x, t)
=D −
∂x x=0 ∂x x=−L
 
∂C(x, t)
=D . (3.3)
∂x x=0

Donde hemos considerado que el flujo en los extremos del cilindro original
es cero, con lo cual el segundo término en la segunda lı́nea del desarrollo
anterior es nulo.

8
 La generalización a (3.3) es la llamada primera ley de Fick :

J~ = −D∇C. (3.4)

A la ecuación de difusión obtenida en la sección anterior [ec. (2.18)], se

le conoce como segunda ley de Fick. Interpretemos el resultado (3.4). Dicha
ecuación nos dice que el flujo de partı́culas ocurre de las zonas con mayor
concentración a las de menor concentración.
La ecuación de continuidad (conservación del número de partı́culas), se
obtiene de forma directa al usar conjuntamente las ecuaciones (2.18) y (3.4).
En efecto, dado que ∇2 = ∇ · ∇, tomando la divergencia en la ec. (3.4) y
usando posteriormente la ec. (2.18), uno encuentra la relación,

∂C(~r, t)
∇ · J~ = − . (3.5)
∂t

Por último hablaremos de la corriente I, ~ concepto que resultará útil en

ejemplos de las subsecuentes secciones. Esta cantidad se define como I~ = AJ. ~
Un análisis rápido de las dimensiones, nos permite establecer que la corriente
mide el número de partı́culas que pasan por una sección plana por unidad de
tiempo. Con la ayuda de la ec. (3.4), podemos expresar la corriente como,

I~ = −AD∇C. (3.6)

Nótese que una vez hallada la concentración, el flujo y la corriente se

obtienen directamente derivando.

4. Formulación de Langevin
En 1908 el fı́sico francés Paul Langevin (1872-1946) presentó una descrip-
ción macroscópica del movimiento browniano. Las consideraciones que hizo
son las siguientes: si en un fluido introducimos una partı́cula grande (en com-
paración con las dimensiones de los constituyentes del fluido), de acuerdo con
la hidrodinámica, esta partı́cula sentirá una fuerza opuesta a su movimiento
y proporcional a su velocidad. A esta fuerza se le conoce como fuerza de fric-
ción viscosa. Por otra parte, las colisiones que experimenta la partı́cula con
las moléculas del fluido, son tomadas en cuenta al introducir una fuerza que
varı́a de forma muy azarosa en el tiempo. Si m representa la masa de esta

9
partı́cula y v su velocidad, la ecuación de movimiento de la misma (conocida
como ecuación de Langevin) es la siguiente,

mv̇ = f (t) − ξv, (4.1)

donde ξ es el coeficiente de fricción y contiene información sobre la viscosi-

dad del fluido, ası́ como sobre la forma geométrica de la partı́cula. La fuerza
estocástica es f (t), y hemos empleado la notación de Newton para la dife-
renciación.
Haremos un breve paréntesis aquı́ para mencionar el significado de una
función de autocorrelación y comentar algo más acerca de la fuerza estocásti-
ca f (t). La función de autocorrelación se define como,
Z ∞
hx(0)x(t)i = dt0 x∗ (t0 )x(t0 + t), (4.2)
−∞

donde el asterisco representa la conjugación compleja. Esta función mide

qué tanta relación existe entre la cantidad x al tiempo t con ella misma
a un tiempo distinto. Ahora, por lo dicho sobre la naturaleza de la fuerza
estocástica f (t), a saber, que fluctúa rápidamente en el tiempo, no esparamos
que su valor a un tiempo t tenga correlación alguna con su valor a otro tiempo.
La manera matemática de expresar esto es la siguiente,

hf (t0 )f (t)i = A δ(t − t0 ), (4.3)

con A una constante con unidades de fuerza al cuadrado a ser determinada.

Nótese también que por su misma naturaleza aleatoria, hf (t)i = 0, para todo
tiempo.
Volvamos a la tarea de encontrar la solución de la ecuación de Lange-
vin. Primero, nótese que en ausencia de la fuerza estocástica, la ec. (4.1) se
convierte en una ecuación diferencial homogénea muy sencilla, cuya solución
está dada por,
v(t) = v(0)e−ξt/m . (4.4)
El tiempo de relajación es m/ξ. Este es el tiempo que le toma a la partı́cula
“olvidar” su posición inicial. Mientras más viscoso es el fluido, más pequeño
será el tiempo con que relaja la velocidad de la partı́cula. Ahora considera-
remos la ec. (4.1) con el término estocástico. Multiplicamos ambos lados de
(4.1) por x, y promediamos sobre un número grande partı́culas,

hmxẍi = hxf (t)i − hξxẋi. (4.5)

10
 Ahora, como
d d
xẋ = xẍ + ẋ2 → xẍ = xẋ − ẋ2 , (4.6)
dt dt
se tiene que,
d
m hxẋi = mhẋ2 i + hxf (t)i − ξhxẋi. (4.7)
dt
Las dos ideas que emplearemos a continuación para seguir con el desarrollo
son: uno no esperarı́a que la posición de la partı́cula esté correlacionada con
la fuerza estocástica en momento alguno, i. e. hxf (t)i = 0, para todo tiempo
t. Por otro lado, el teorema de equpartición asienta que,
1 1
mẋ2 = kB T, (4.8)
2 2
donde kB es la constante de Boltzmann. Usando lo anterior, la ec. (4.7) queda
pues como sigue,
d
m hxẋi = kB T − ξhxẋi. (4.9)
dt
Vemos que la ecuación resultante es una ecuación diferencial lineal no ho-
mogénea de primer orden. Es fácil comprobar que la solución a esta ecuación
está dada por,
kB T
hxẋi = C0 e−ξt/m + , (4.10)
ξ
donde C0 se determina imponiendo una condición inicial. En nuestro caso,
seguiremos considerando que las partı́culas inician su movimiento desde el
origen, con lo cual hxẋi = 0 en t = 0, de donde resulta que C0 = − kBξ T , y
llegamos a que,
kB T
hxẋi = (1 − e−ξt/m ). (4.11)
ξ
Luego, usando la identidad
1d 2
x = xẋ, (4.12)
2 dt
nos resulta la ecuación,
1d 2 kB T
hx i = (1 − e−ξt/m ), (4.13)
2 dt ξ

11
 Analicemos ahora el comportamiento asintótico de la solución para tiem-
pos muy pequeños y muy grandes, respectivamente, comparados con el tiem-
po de relajación m/ξ, asociado a la velocidad. Si t  m/ξ, desarrollamos en
serie de Taylor la exponencial en la expresión anterior,
ξ
e−ξt/m ≈ 1 − t,
m
con lo cual
d 2 2kB T ξ 2kB T
hx i = (1 − (1 − t)) = t. (4.14)
dt ξ m m
Al integrar esta expresión, obtenemos que en este lı́mite,
kB T 2
hx2 i = t. (4.15)
m
Si ahora t  m/ξ, la exponencial en (4.13) es prácticamente nula, y al
integrar la expresión resultante, obtenemos,
2kB T
hx2 i = t. (4.16)
ξ
Como se puede observar, esto corresponde a una partı́cula que difunde
siguiendo una caminata aleatoria simple (ver sección 2), donde identificamos
el coeficiente de difusión,
kB T
D= . (4.17)
ξ
Nótese que la constante de difusión D no depende de la masa de la partı́cu-
la browniana.
Antes de avanzar, vale la pena esclarecer una cuestión interesante que sur-
ge del desarrollo antes hecho. Si repasamos con cuidado las últimas lı́neas,
notaremos primero que hemos utilizado un resultado que invoca a una situa-
ción de equilibrio térmico. En efecto, para llegar a la ec. (4.9) empleamos el
teorema de equipartición. Uno podrı́a entonces pensar que el resultado al que
se llega considerando tiempos pequeños [ec. (4.15)] es inconsistente con esta
manera de proceder, pues si se trata de una situación en que inclusive las ve-
locidades no han relajado, obviamente se está lejos de alcanzar el equilibrio.
Sin embargo, el lı́mite (4.15) es correcto y de hecho tiene bastante sentido
fı́sico. Veamos porqué esto es ası́.

12
 La clave está en identificar con precisión las diferentes escalas de tiempo
en el comportamiento de la partı́cula browniana. Primero, para tiempos mu-
cho menores que m/ξ, el tiempo de relajación de la velocidad, la partı́cula
browniana deberı́a comportarse como partı́cula libre, pues hablamos de tiem-
pos tan cortos que aún no siente que está inmersa en el fluido. En tal caso, la
raı́z cuadrada del desplazamiento cuadrático medio debe ser proporcional a
t (velocidad constante), tal como lo indica la ec. (4.15). Luego, para tiempos
mucho mayores que m/ξ pero no tan “grandes”(y en seguida aclararemos
a qué nos referimos con esta expresión), la partı́cula entra al llamado régi-
men difusivo. La interacción con el fluido hace que el comportamiento de la
partı́cula se modifique completamente, y ahora el desplazamiento de ésta va
de acuerdo con la ec. (4.16). Aquı́ se llega a una situación de equilibrio con
el fluido, al menos en lo que a velocidades respecta. Por último tenemos el
régimen de tiempos hidrodinámicos. Para estos tiempos, la manera en que
las partı́culas se reparten en el espacio es homogénea, de forma que se tiene
ahora un equilibrio total. Estos tiempos, para casos fı́sicos reales, suelen ser
mucho mayores que el tiempo de relajación de la velocidad (la comparación
involucra nueve órdenes de magnitud). Tiempos menores que el tiempo hi-
drodinámico, son los tiempos no tan “grandes”, a los que nos referimos en el
régimen difusivo.
Entonces, que se utilice el teorema de equipartición cuando t  m/ξ,
está bien justificado. Pero lo que podemos concluir ahora es que para t 
m/ξ, estamos en otra situación en la que es posible emplear el mismo resul-
tado. En efecto, el movimiento traslacional de un gas monoatómico, es quizás
el ejemplo más sencillo para ilustrar el teorema de equipartición, y es dicha
situación la que se tiene en la descripción del movimiento browniano para
tiempos pequeños. De manera que, de haber hecho la distinción clara entre
los comportamientos tan distintos de la partı́cula browniana dependiendo de
la escala de tiempo considerada, se justifica el empleo del teorema de equi-
partición en los lı́mites que llevan a las ecuaciones (4.15) y (4.16), por lo que
ambas expresiones son correctas en los dominios correspondientes.
Por último, calcularemos la función de correlación de la velocidad, y en-
contraremos la relación del coeficiente de fricción ξ, con la función de corre-
lación de la fuerza f (t). Para calcular la correlación de la velocidad, multi-
plicamos por v(0) la ec. (4.1) y el promediamos,
D d E
m v(0) v(t) = hv(0)f (t)i − ξhv(0)v(t)i. (4.18)
dt

13
 Nótese que el primer término del lado derecho en la ecuación anterior
es nulo. En efecto, hv(0)f (t)i = v(0)hf (t)i = 0. Con lo cual, nos queda
una ecuación diferencial muy sencilla para la función de correlación de la
velocidad,
d ξ
hv(0)v(t)i + hv(0)v(t)i = 0, (4.19)
dt m
cuya solución viene dada por,

hv(0)v(t)i = hv(0)2 ie−ξt/m . (4.20)

Esta ecuación nos dice que mientras más tiempo pase, menor será la
correlación entre velocidades.
Finalmente, relacionaremos la constante A que aparece en (4.3) con la
temperatura T . Para esto, comenzamos reescribiendo la ecuación de Langevin
(4.1), de la siguiente manera,

d ξ f (t)
v(t) + v(t) = ≡ F (t). (4.21)
dt m m
La ecuación (4.21), es una ecuación diferencial lineal no-homogénea de
primer orden. La solución a este tipo de ecuaciones puede ser encontrada
usando el método del factor integrante. La idea básica de este método, con-
siste en multiplicar la ecuación diferencial por una función µ(t) tal que el
lado izquierdo de (4.21) pueda escribirse como la derivada con respecto a t
del producto µ(t)v(t), para después integrar directamente la expresión resul-
tante. Se comprueba fácilmente que si µ(t) = exp(ξt/m), es posible hacer lo
planteado lı́neas arriba. Entonces, (4.21) toma la forma,

d ξt/m
(e v(t)) = eξt/m F (t). (4.22)
dt
Integrando esta ecuación, de t = 0 hasta un tiempo arbitrario t, uno
encuentra que, Z t
ξt/m
e v(t) − v(0) = dτ eξτ /m F (τ ),
0
de donde se sigue que
Z t
−ξt/m
v(t) = v(0)e + dτ eξ(τ −t)/m F (τ ). (4.23)
0

14
 Teniendo en mente la ec.(4.3) y el hecho de que hf (t)i = 0, elevamos
al cuadrado ambos miembros de la ec.(4.23) y promediamos posteriormente.
Encontraremos un término cruzado en el lado derecho que involucra al pro-
medio hF (t)i, por lo que dicho término es nulo. Otro de los términos contiene
el producto,
Z t Z t
−ξ(t−τ )/m 0
dτ e F (τ ) dτ 0 e−ξ(t−τ )/m F (τ 0 ), (4.24)
0 0

al promediar este término, por (4.3) resulta lo siguiente,

Z t Z t
−ξ(t−τ )/m 0 A
dτ e dτ 0 e−ξ(t−τ )/m 2 δ(τ − τ 0 ) =
0 0 m
Z t
A A
= 2 dτ e−2ξ(t−τ )/m = (1 − e−2ξt/m ). (4.25)
m 0 2mξ
Finalmente, tenemos que,
A
hv 2 i = v(0)2 e−2ξt/m + (1 − e−2ξt/m ). (4.26)
2mξ
Para tiempos muy largos, las dos exponenciales en lado derecho de la
ecuación anterior son despreciables, y si dichos tiempos son tales que se ha
llegado ya al equilibrio térmico, el teorema de equipartición vale, con lo que
resulta la relación,
kB T A
= , (4.27)
m 2mξ
de donde se sigue que,
A = 2ξkB T. (4.28)
Este resultado es una forma del teorema de fluctuación-disipación. Rela-
ciona la magnitud A, de la fuerza fluctuante f (t), con la fricción o disipación
contenida en el coeficiente ξ. Utilizando este resultado en (4.3), se sigue que,
Z ∞
1
ξ= dt hf (0)f (t)i. (4.29)
2kB T −∞

15
5. Solución general a la ecuación de difusión
En esta sección estudiaremos uno de los métodos clásicos de la fı́sica
matemática para encontrar la solución general a la ecuación de difusión. El
método en cuestión es el método de la transformada de Fourier. La idea detrás
de esta técnica, consiste en pasar de una ecuación en derivadas parciales a
una ecuación diferencial ordinaria dependiente sólo del tiempo.
Dada una función f (~r), definimos su transformada de Fourier como,
Z ∞
f˜(ω) = F{f (~r)} := d3~r f (~r)e−iω~r . (5.1)
−∞

Para lo que sigue, haremos uso de dos propiedades de la transformada de

Fourier que se prueban fácilmente a partir de la definición (5.1). La primera
propiedad es la de linealidad,

F{f (x) + g(x)} = F{f (x)} + F{g(x)} = f˜(x) + g̃(x), (5.2)

y la segunda, es la transformada de la derivada, la cual se prueba integrando

por partes,
n df o
F = iω f˜(x). (5.3)
dx
Con estas dos propiedades, es sencillo establecer que el resultado de apli-
car la transformada de Fourier a la ecuación de difusión con arrastre (por
simplicidad, consideraremos de momento el caso unidimensional), es el si-
guiente,
∂ C̃(ω, t)
= −iων C̃(ω, t) − ω 2 DC̃(ω, t)
∂t
= −(iων + ω 2 D)C̃(ω, t). (5.4)

Esta es una ecuación diferencial de variables separables y se puede integrar

fácilmente. La solución a (5.4), está dada por,
2 D)t
C̃(ω, t) = C̃0 e−(iων+ω . (5.5)

Para expresar la solución en términos de la variable x, aplicamos la trans-

formada inversa de Fourier a la ec. (5.5), la cual se define como,
Z ∞
1
f (~r) = dω f˜(ω)eiω~r . (5.6)
2π −∞

16
 Se obtiene que,
Z ∞
1 2 D)t
C(~r, t) = dω C̃0 e−(iων+ω eiω~r . (5.7)
2π −∞

La integral del lado izquierdo aparece muy frecuentemente en los libros

de texto de variable compleja, como ejemplo de una aplicación del teorema
de Cauchy para evaluar integrales.4 Para no extender en deması́a el análisis
presentado en esta sección, daremos únicamente el resultado al que se llega
después de evaluar dicha integral.
C0 n (x − νt)2 o
C(x, t) = √ exp − . (5.8)
4πDt 4Dt
Esta distribución es gaussiana 5 con desviación cuadrática media σ 2 =
2Dt. La solución de la ecuación de difusión sin arrastre se obtiene poniendo
ν = 0 en la expresión anterior. Se obtiene ası́, una solución al problema de
difusión en toda la recta real, i.e. para −∞ < x < ∞.
Una de las propiedades de invarianza básicas de la ecuación de difusión
es que una integral de soluciones de la misma, es de hecho solución. Esto
es, si G(x, t) es una solución a la ecuación de difusión, entonces G(x − x0 , t)
también lo es (ya que al hacer una traslación del argumento no se modifica
el carácter de una solución), y lo es también la integral,
Z ∞
dx0 G(x − x0 , t)g(x0 ), (5.9)
−∞
0
para alguna función g(x ) tal que esta integral impropia converga apropiada-
mente. Se encuentra que (5.9) da la solución única al problema de difusión,
para una distribución inicial dada por la función g(x), siempre que,6
1 n x2 o
G(x, t) = √ exp − , (5.10)
4πDt 4Dt
4
Ver por ejemplo: Stein EM, Shakarchi R, Complex Analysis (Princeton Lectures in
Analysis, No. 2), Nueva Jersey: Princeton University Press, 2003, pp. 42–44.
5
Una distribución gaussiana con media µ y desviación estándar σ, viene dada por:
1 n (x − µ)2 o
N [µ, σ] = √ exp − .
σ 2π 2σ 2
6
Los detalles matemáticos detrás de esta afirmación no son tan complicados y pueden
ser consultados en: Strauss WA, Partial differential equations: an introduction, Estados
Unidos: John Wiley & Sons, Inc., 1992, pp. 45–50.

17
para todo t > 0. La función G(x, t) se conoce como función de Green. En
definitiva, conocida la concentración inicial C(x, 0) = C0 , la solución general
al problema de difusión sobre toda la recta real, viene dada por,
Z ∞ n (x − x0 )2 o
C0
C(x, t) = √ dx0 exp − . (5.11)
4πDt −∞ 4Dt

6. Condiciones de frontera
El siguiente paso en nuestro estudio de la ecuación de difusión es, natu-
ralmente, analizar las situaciones donde se imponen condiciones de contorno
sobre la concentración para valores finitos de las coordenadas espaciales.
Cuando tratemos un problema en el que sea necesario satisfacer semejan-
tes condiciones, tendremos de hecho lo que se llama, un sistema confinado.
Particularmente, nos centraremos en el estudio de cuatro tipos de fronteras, a
saber: fronteras absorbentes, reflejantes, parcialmente absorbentes, o cuando
éstas aportan partı́culas al sistema.
Una pared absorbente remueve del sistema cualquier partı́cula que entre
en contacto con ella. La concentración de partáculas sobre toda la pared, es
cero para todo tiempo. Si ~r es un punto de la misma, entonces,
C(~r, t) = 0, para todo t > 0, (6.1)
da la condición matemática para representar un pared de este tipo.
Se dice que la frontera es reflejante, cuando las partı́culas cambian de
dirección al entrar en contacto con ésta. El flujo de partı́culas a través de la
misma es cero, y esto se representa de la siguiente manera,
∇C(~r, t) = 0, para todo t > 0. (6.2)
Luego tenemos un caso intermedio entre los dos anteriores: las paredes
parcialmente absorbentes. Cuando se quiera modelar este tipo de frontera,
uno introduce una constante κ que da cuenta de la eficiencia con que la
frontera permite el paso de partı́culas. Esta condición se expresa como sigue,
∇C(~r, t) = κC(~r, t) para todo t > 0. (6.3)
El último caso del que nos ocuparemos, es el de paredes que aportan
partı́culas al sistema. Esto puede ocurrir si la naturaleza quı́mica de las fron-
teras es tal que tienen lugar reacciones que son las que suministrarán las

18
nuevas partı́culas. Si la susodicha reacción se lleva a cabo con una constante
de reación k, entonces, la condición que se debera imponer es,
∂C(~r, t)
= kC(~r, t). (6.4)
∂t
La solución de (6.4) da la forma correcta de la concentración sobre la
forntera para todo tiempo. Luego, se resuelve el la ecuación de difusión im-
poniendo una condición de frontera tipo del Dirichlet no homogénea.

7. Solución a la ecuación de difusión: Estado

estacionario
En esta sección estudiaremos algunos ejemplos de difusión en geometrı́as
sencillas. Lo haremos además, bajo el supuesto de que el sistema se encuentra
en un estado estacionario, lo cual quiere decir que la concentración no cambia
más en el tiempo, esto es, ∂C/∂t = 0. Con lo cual, la ecuación de difusión se
reduce a la ecuación de Laplace:
∇2 C(~r, t) = 0. (7.1)
Si recordamos la ecuación de continuidad [ec. (3.5)], vemos también que
en el régimen estacionario ocurre que ∇ · J~ = 0. Cuando además de estar en
un estado estacionario el flujo es nulo, J~ = ~0, se define el estado de equilibrio,
pues en tal caso por (3.4) tenemos que ∇C = 0, lo cual implica que la
concentración permanece constante en el espacio y tiempo.

7.1. Difusión a través de una membrana

La geometrı́a del problema que vamos a tratar se muestra en la Fig. 1.
Se trata pues de una membrana de espesor h con un extremos en x = 0, h.
Supondremos que las concentraciones en los contenedores separados por esta
membrana son Ci y Cd , en el lado izquierdo y derecho, respectivamente, y
que Ci > Cd , de manera que flujo va de izquierda a derecha (primera ley de
Fick). Ahora bien, dada la simetrı́a radial del problema podemos asumir que
la concentración sólo depende de la coordenada x, con lo cual, al suponer un
flujo estacionario, la ecuación que determina el perfil de concentración es:
d2 C(x)
= 0. (7.2)
dx2
19
Figura 1: Difusión a través de una membrana. Suponemos que la concentración Ci
del compartimento en el lado izquierdo de la imagen es mayor que la concentración
Cd del volumen derecho, razón por la cual el flujo se da en la dirección indicada
por las flechas. Colocamos el sistema de manera que el eje del cilindro coincida
con el eje x.

Integrando consecutivamente dos veces, obtenemos que,

C(x) = A1 x + A2 , (7.3)

donde A1 , A2 , son constantes de integración a ser determinadas. La ec. (7.3)

muestra que el perfil de concentración es lineal en x en el interior de la
membrana. Determinaremos ahora las constantes de integración. Para esto,
supondremos que la membrana permite el paso de partı́culas en cualquiera
de sus dos caras con una eficacia cuantificada con la constante K. Ası́, la
concentración en x = 0 es,

C(x = 0) = KCi . (7.4)

De igual manera la concentración en la cara derecha de la membrana es,

C(x = h) = KCd . (7.5)

Imponiendo estas condiciones a (7.3), uno encuentra que A2 = KCi y que

A1 = K(Cd − Ci )/h. Esto es,

(Cd − Ci )
C(x) = K x + KCi . (7.6)
h
El comportamiento de esta cantidad al variar K manteniendo h constante
y viceversa, está ilustrado en la Fig. 2.El flujo de partı́culas a través de la
membrana puede calcularse usando la ec. (3.4):
dC DK
J = −D = (Ci − Cd ). (7.7)
dx h
20
 Como uno podrı́a esperar intuitivamente, el flujo aumenta de forma pro-
porcional con K y disminuye al incrementar el espesor de la membrana h.
La permeabilidad de la membrana se define como el flujo divido por la di-
ferencia de concentración en ambos lados de la membrana. En este caso, la
permeabilidad es DK/h.
HaL cd =0.01ci , h=3e-8m HbL cd =0.01ci , K=3e-8m
4. ´ 10-8

8. ´ 10-8
Kh=1 Kh=1
3. ´ 10-8
-8
6. ´ 10 Kh=1.5 Kh=1.5

Kh=2
cHxLci

Kh=2

cHxLci
2. ´ 10-8
4. ´ 10-8

1. ´ 10-8
2. ´ 10-8

0 0
0 5. ´ 10-9 1. ´ 10-8 1.5 ´ 10-8 2. ´ 10-8 2.5 ´ 10-8 3. ´ 10-8 0 1. ´ 10-8 2. ´ 10-8 3. ´ 10-8 4. ´ 10-8 5. ´ 10-8 6. ´ 10-8
x x

Figura 2: Gráficas del perfil de concentración de partı́culas difundiendo a través

de una membrana. (a) Se muestra el efecto mantener el ancho de la membrana h
constante, variando K. (b) Se muestra el efecto de variar h, para una K dada.

7.2. Difusión hacia una proteı́na o micela

En muchos problemas de la biologı́a celular interesa estudiar el proceso de
captura por difusión, i.e. el proceso mediante el cual especies que difunden
llegan a cierto cuerpo, el cual las absorbe. Tal puede puede ser el caso de
ligandos que son capturados por receptores distribuidos sobre una superficie
esférica. Una primera aproximación para abaordar el problema anterior, es
suponer que en el medio donde se encuentran los ligandos introducimos una
esfera con propiedades perfectamente absorbentes en su frontera (ver Fig. 3).
Por la simetrı́a del problema, estableceremos como hipótesis adicional que
la concentración sólo tiene dependencia radial. Con esta supoción y usando
coordenadas esfeéricas,7 la ecuación de Laplace (seguimos considerando es-
7
En coordenadas esféricas el operador laplaciano es:

1 ∂  2 ∂  1 ∂  ∂  1 ∂2
∇2 = r + senθ + .
r2 ∂r ∂r r2 senθ ∂θ ∂θ r2 sen2 θ ∂φ2

21
 Figura 3: Difusión hacia una esféra con superficie de naturaleza absorbente.

tados estacionarios) adopta la forma siguiente:

1 d2 (rC)
∇2 C = = 0, (7.8)
r dr2
la cual implica que,
d(rC)
= A1 , (7.9)
dr
con A1 una constante de integración. Integrando nuevamente obtenemos,

rC(r) = A1 r + A2 ,

o bien,
A2
C(r) = A1 + . (7.10)
r
Sigue ahora determinar las constantes de integración. Para ello, además
de la condición sobre la superficie de la esfera:

C(r = a) = 0, (7.11)

donde a es el radio de la misma, necesitamos una condición adicional. Esta

condición tiene que ver con el valor de la concentración muy lejos de la esfera.

22
Impondremos que en tales puntos tenga un valor constante C∞ , esto es,
C(r → ∞) = C∞ . (7.12)
Con estas restricciones, es inmediato ver que A1 = C∞ y A2 = −aC∞ ,
con lo cual obtenemos,  a
C(r) = C∞ 1 − . (7.13)
r
El comportamiento de esta solución se muestra en la Fig. 4 para distintos
valores de a. El flujo de partı́culas está dado por,
dC aDC∞
J = −D =− . (7.14)
dr r2

1.0

0.8

0.6
CHrLC¥

a=1
0.4

a=1.5

0.2 a=2

0.0
0 1 2 3 4 5 6 7 8 9 10
r
a

Figura 4: Gráficas del perfil de concentración de partı́culas difundiendo hacia una

esfera de radio a y con superficie de naturaleza absorbente. Se muestra el efecto
de variar el radio de la esfera.

Multiplicando este flujo por la superficie de la esfera 4πa2 , uno tiene la

corriente, la cual da el número de colisiones que tienen lugar en la superficie
de la esfera por unidad de tiempo.
I = −4πDC∞ a. (7.15)

23
7.3. Dfusión acoplada a una reacción quı́mica
El último ejemplo que veremos en esta sección es el de partı́culas capaces
de reaccionar quı́micamente con el medio en el cual difunden. El hecho de
que las partı́culas puedan transformarse quı́micamente, implica que la con-
centración disminuirá al ritmo dictado por la constante de reacción kr (en
unidades de inverso de tiempo). Al incluir un proceso de reacción, la ecuación
que gobierna la concentración es ahora la siguiente:

∂C ∂ 2C
= D 2 − kr C, (7.16)
∂t ∂x
donde D es la constante de difusión del medio. Considerando el caso estacio-
nario, la ecuación anterior se reduce a,

∂ 2 C kr
− C = 0, (7.17)
∂x2 D
que es una ecuación diferencial ordinaria de segundo orden elemental. La
solución general a (7.17) viene dada por,

C(x) = A1 eαx + A2 e−αx , (7.18)

p
con α := kr /D, y A1 , A2 , constantes cuyo valor se elije de manera que
de cumplan las condiciones de frontera del problema. Las condicones que
aquı́ tomaremos son las siguientes: (i) a muy grandes distancias del origen,
la concentración será nula, C(x → ∞) = 0. (ii) En el origen del sistema de
coordenas, establecemos que se tenga una concentración constante C0 . De
imponer la primera de estas condiciones, resulta A1 = 0. Al hacer lo propio
con la segunda, tenemos que A2 = C0 . Con lo cual la solución al problema
en cuestión queda dada por:

C(x) = C0 e−αx . (7.19)

La concentración decae exponencialmente con la distancia, haciéndolo

más rápido entre mayor sea la constante de reacción kr y menor el coeficiente
de difusión D. Este comportamiento se muestra en la Fig. 5 El flujo en el
origen es, p
J|x=0 = DC0 α = C0 kr D. (7.20)

24
 1.0

Α=1
0.8

Α=0.5

0.6
Α=0.2
CHxLC0

0.4

0.2

0.0
0 2 4 6 8 10
x

Figura 5: Concentración de partı́culas

p en el caso de difusión acoplada a una reacción
quı́mica. El parámetro α = kr /D, mide el efecto conjunto de la tasa de reacción
y de la difusión.

8. Probabilidad de sobrevivencia y tiempo

promedio de captura
Dos conceptos que resultan de suma utilidad en el análisis de la difusión
en sistemas confinados son la probabilidad de sobrevivencia y el tiempo pro-
medio de captura. Antes de proseguir, resaltaremos una simetrı́a a la que no
hemos hecho alusión antes. El problema de calcular la probabilidad de en-
contrar una partı́cula browniana, digamos en el punto x al tiempo t, siendo
que la partı́cula se hallaba inicialmente en x0 , es análogo al de preguntar la
probabilidad de que dada una partı́cula en x al tiempo t, ésta haya iniciado
su caminata en el punto x0 . Esta simetrı́a, consecuencia de que el movimien-
to browniano sigue la propiedad de Markov, implica que si cambiamos en
la ecuación de difusión las derivadas parciales con respecto a x por las co-
rrespondientes con respecto a x0 , la probabilidad condicional P (x, t|x0 ), de
hallar a la partı́cula en x al tiempo t, dado que inició en x0 , también ha de

25
satisfacer la ecuación resultante, esto es,
∂P (x, t|x0 ) ∂ 2 P (x, t|x0 )
=D (8.1)
∂t ∂x20
Por lo descrito lı́neas arriba, a esta ecuación se le conoce como ecuación
hacia atrás, y nos servirá para calcular las propiedades antes mencionadas.
Comenzaremos por considerar el caso en que las fronteras que delimitan
nuestro sistema son perfectamente absorbentes. Nótese que tiene mucho sen-
tido estudiar la probabilidad de sobrevivencia y el tiempo medio de captura
con la ec. (8.1) [i.e. considerando variaciones con respecto a la posición ini-
cial x0 ], pues claramente disminuirán sus valores mientras más cerca de la
frontera se coloque la partı́cula inicialmente.
La probabilidad de sobrevivencia se define como,
Z L
S(t|x0 ) := dx P (x, t|x0 ), (8.2)
0

donde L es una longitud que caracteriza la dimensión del sistema. Como

puede observarse, esta cantidad nos da la probabilidad de encontrar a la
partı́cula en cualquier punto del sistema al tiempo t. Para determinar la
ecuación que la gobierna, integramos de 0 a L la ec. (8.1). Haciendo esto uno
obtiene, para el lado izquierdo,
Z L
∂ L
Z
∂P (x, t|x0 ) ∂S(t|x0 )
dx = dx P (x, t|x0 ) = ,
0 ∂t ∂t 0 ∂t
y para el miembro derecho,
Z L Z L
∂ 2 P (x, t|x0 ) ∂2 ∂ 2 S(t|x0 )
dx D = D dx P (x, t|x 0 ) = D .
0 ∂x20 ∂x20 0 ∂x20
Esto es, la probabilidad de sobrevivencia S(t|x0 ), es gobernada por la
misma ecuación que rige a P (x, t|x0 ):
∂S(t|x0 ) ∂ 2 S(t|x0 )
=D . (8.3)
∂t ∂x20
Una vez definida la probabilidad de sobrevivencia, se define el tiempo
promedio de captura como,
Z ∞
τ (x0 ) := dt S(t|x0 ), (8.4)
0

26
 Esta cantidad nos da la probabilidad de que la partı́cula sobreviva para
todo tiempo. Para encontrar la ecuación la gobierna, procedemos como antes
e integramos con respecto a t la ec. (8.3) de 0 a ∞. Esto da,
Z ∞ Z ∞
∂S(t|x0 ) ∂ 2 S(t|x0 )
dt = dt D ,
0 ∂t 0 ∂x20
Z ∞ Z ∞
∂2
dS(t|x0 ) = D 2 dt S(t|x0 ),
0 ∂x 0
esto es,
∂ 2τ
S(t|x0 )|∞
0 = D . (8.5)
∂x20
Si inicialmente la partı́cula se encuentra en el seno del sistema, entonces
S(t = 0|x0 ) = 1. Luego, la partı́cula eventualmente alcanzará la frontera, con
lo cual, para tiempos muy largos, S(t → ∞|x0 ) = 0. Sustuituyendo esto en
(8.5) encontramos la ecuación que satisface τ (x0 ), a saber,

∂ 2 τ (x0 ) 1
2
=− . (8.6)
∂x0 D

En tres dimensiones (8.6) adopta la forma de la ecuación de Poisson:

1
∇2 τ (~r0 ) = − , (8.7)
D

con ~r0 el vector que da la posición inicial de partı́cula.

A continuación daremos un ejemplo donde veremos como aplicar las ideas
recién desarrolladas. Considérese una partı́cula que difunde en una dimen-
sión. La longuitud de la celda en la que difunde es L, y sus paredes se en-
cuentran colocadas en x = 0, L, siendo ambas de naturaleza absorbente. Para
calcular el tiempo promedio de sobrevivencia de bbemos resolver la ec. (8.6)
sujeta a las condiciones mencionadas, a saber,

τ (x0 = 0) = 0; τ (x0 = L) = 0. (8.8)

La solución general de (8.6) es,

1 2
τ (x0 ) = − x + A1 x0 + A2 , (8.9)
2D 0

27
que al imponer las condiciones (8.8) deriva en,
1 2 L x0
τ (x0 ) = − x0 + x0 = (L − x0 ). (8.10)
2D 2D 2D
Recordemos que τ da el tiempo que en promedio le toma a la partı́cula
escapar del sistema, cuya longitud es L. Podemos hacer dos variaciones del
ejercicio anterior. Primero, supongamos que la partı́cula puede iniciar en
cualquier punto del intervalo [0, L] con la misma probabilidad. En tal caso,
debemos promediar (8.10) sobre toda la recta de longitud L, esto es,
1 L x0 (L − x0 )
Z
hτ (x0 )i = dx0
L 0 2D
3 3i
1 Lh L L2
= − = . (8.11)
2DL 2 3 12D
La segunda variación que podemos estudiar de forma inmediata, es pensar
que una de las paredes, digamos la que se encuentra en x0 = L, es reflejante,
con lo cual en ese extremo la condición a ser impuesta es ahora,
∂τ
= 0. (8.12)
∂x0 x0 =L

Las constantes de integración que aperecen en (8.9), resultan ser en este

caso A1 = L/D y A2 = 0, con lo que,
x0  x0 
τ (x0 ) = L− . (8.13)
D 2
Por último, vale la pena resaltar que las ecs. (8.6, 8.7), nos brindan un
medio relativamente sencillo para estimar coeficientes de difusión de forma
empı́rica para sistemas con diferentes condiciones iniciales y de frontenra,
teniendo que calcular únicamente los tiempos promedios de sobrevivencia.

9. Procesos de difusión bajo confinamiento:

Sistemas formados por dos cavidades in-
terconectadas
Resolver la ecuación de difusión [ec. (2.18)] con condiciones iniciales y a
la frontera arbitrarias es una tarea muy difı́cil de cumplir e inclusive algunas

28
veces resulta prácticamente imposible. En general, no será una actividad
trivial obtener la función C(~r, t), y debemos hacernos de metodologı́as que
nos permitan extraer información del sistema bajo estudio. Uno de estos
enfoques se presentó en la sección anterior, donde se mostró cómo calcular
el tiempo promedio del primer arribo a un cierto dominio. Tal estrategia
sirve para elucidar al menos la parte dinámica del proceso de difusión en
geometrı́as complejas.
En esta sección estudiaremos el problema de difusión en sistemas forma-
dos por dos cavidades interconectadas por un orificio de forma circular de
radio a. No pretendemos encontrar la concentración C(~r, t), para cualquier
punto y en todo momento, sino los cambios de esta distribución en los dos
volúmenes considerados.
Estudiar cómo difunde una partácula browniana entre dos repertorios es
análogo a revisar un problema de un sistema capaz de ocupar únicamente
dos estados y analizar las transiciones entre éstos.
La expresión que relaciona la corriente que pasa a través de un disco con
el radio del mismo es conocida y se debe a Terrell L. Hill. Dicha fórmula es
la siguiente,

4aD
k= , (9.1)
V
donde a es el radio del disco, D es la constante de difusión y V el volumen
del contenedor de donde proviene la partı́cula.
La cantidad representada en (9.1) como ha sido dicho es la corriente que
pasa por el disco, o sea, el número de partı́culas que llegan al disco por
unidad de tiempo y puede asociarse con una constante de velocidad en una
cinética de primer orden (este término está referido a una reacción donde una
sustancia A se transforma en otra sustancia B sin la formación de productos
intermediarios). Veamos cómo se deduce el resultado (9.1).

9.1. Fórmula de Hill

Tenemos el problema de la difusión de partı́culas hacia un disco absor-
bente de radio a embebido en un plano reflejante, digamos el X–Y . La ecua-
ción que describe este fenómeno como hemos visto es la ecuación de Fick
[ec. (2.18)]. Consideraremos además el caso estacionario, donde la concentra-
ción es independiente del tiempo, de manera tal que la ecuación que realmente

29
debemos resolver es la ecuación de Laplace:

∇2 C = 0. (9.2)

Dada la simetrı́a del problema, es conveniente el uso de coordenadas

cilı́ndricas, con lo que la ecuación anterior adopta la forma,
 2
1 ∂2 ∂2

∂ 1 ∂
+ + + C = 0. (9.3)
∂r2 r ∂r r2 ∂θ2 ∂z 2

Nuevamente, teniendo en mente a la simetrı́a de la situación considerada,

supondremos que la concentración sólo depende de r y z. Luego, buscamos
una solución de (9.3) de la forma

C(r, z) = R(r)Z(z). (9.4)

Sustituyendo (9.4) en (9.3) y dividiendo ambos lados de la ecuación re-

sultante por C, llegamos a

1 ∂ 2 R 1 ∂R 1 ∂ 2Z
 
+ + = 0. (9.5)
R ∂r2 r ∂r Z ∂z 2

Lo cual implica que las funciones R(r) y Z(z), deben satisfacer separa-
damente las ecuaciones,

1 ∂ 2 R 1 ∂R
 
2
+ = −γ 2 , (9.6)
R ∂r r ∂r
1 ∂ 2Z
= γ 2, (9.7)
Z ∂z 2
donde γ 2 es una constante de separación.
Ahora discutiremos cuáles son las condiciones a la frontera que mejor
modelan al sistema en cuestión. Tales condiciones con sus justificación se
mencionan en seguida.

1. Como el disco que consideramos es perfectaente absorbente, la concen-

tración en éste ha de ser nula, esto es,

C(r, 0)|r≤a = 0. (9.8)

30
 2. El resto del plano donde está contenido el disco lo pensamos como
reflejante, por lo que el flujo en en tal extensión del espacio es cero, o
sea,
∂C
= 0. (9.9)
∂z r>a,z=0
3. Ahora, consideraremos que la concentración en el seno del disolvente
tiene un valor constante, i.e. en regiones muya alejadas del disco ab-
sorbente. Entonces, primero imponemos que la concentración sobre el
plano para valores muy grandes de r sea constante, digamos que tome
el valor C∞ . En términos matemáticos esto quiere decir que,

C(r, z)|r→∞,z=0 = C∞ . (9.10)

4. Ası́ mismo, con z muy grande la concentración ha de ser constante para

cualquier valor del radio, es decir,

C(r, z)|r≥0,z→∞ = C∞ . (9.11)

Establecidas las condiciones de contorno, pasemos a resolver las ecuacio-

nes (9.6) y (9.7). La solución a la ec. (9.7) es bien conocida, y está dada
por,
Z(z) = Aeγz + Be−γz . (9.12)
De acuerdo a la última condicón enunciada [ec. (9.11)], la concentración
ha de ser finita para z → ∞, de donde se sigue que A = 0. Con lo cual queda
la solución en la coordenada z:

Z(z) = Be−γz . (9.13)

La determinación de la constante B se hará más adelante. Continuemos

de momento con la solución a la parte radial. La ec. (9.6) puede reescribirse
de la siguiente manera:
∂ 2R ∂R
r2 2
+r + γ 2 r2 R = 0, (9.14)
∂r ∂r
La forma de la expresión resultante, recuerda a un tipo particular de
ecuaciones diferenciales ordinarias lineales, cuya aparición es muy usual en
los problemas de la fı́sica matemática, nos referimos a la ecuación paramétrica
de Bessel :

31
 x2 y 00 (x) + xy 0 (x) + (λ2 x2 − ν 2 )y(x) = 0, (9.15)
que tiene la solución general:

y(x) = c1 Jν (λx) + c2 Yν (λx), (9.16)

con Jν (λx) y Yν (λx), las funciones de Bessel de primera y segunda clase,

respectivamente, y de orden ν.
Identificando γ = λ y ν = 0, tenemos que la solución a la ec. (9.14) es de
la forma,
R(r) = DJ0 (γr) + EY0 (γr). (9.17)
Se conoce que la función Y0 (γr) diverge en el lı́mite r → 0. Como estamos
en la búsqueda de una solución acotada en tal lı́mite, concluimos que la
constante E tiene que ser cero. ası́ nos queda la solución:

R(r) = DJ0 (γr). (9.18)

Sustituyendo (9.13) y (9.18) en (9.4), uno encuentra que,

C(r, z) = F J0 (γr)e−γz , (9.19)

donde F = BD.
Para que las condiciones (9.8) y (9.9) sean satisfechas simultáneamente,
se requiere que la solución sea tal que

J0 (γr) = 0, para r ≥ 0. (9.20)

De manera tal que el parámetro γ sólo puede tomar valores dentro del
conjunto γn de las raı́ces de (9.20). La superposición continua de las soluciones
obtenidas es a su vez solución de (9.3), esto es, la solución formal al problema
planteado está dada por,
Z ∞
C(r, z) = F dγ J0 (γr)e−γz f (γ), (9.21)
0

con f (γ) una función arbitraria del parámetro γ. La forma de esta función,
se escoge teniendo en cuenta la forma de algunas de las integrales definidas
de Bessel. En particular, considerando
sen(γa)
f (γ) = , (9.22)
γ

32
se obtiene,
Z ∞ 
sen(γa) arc sen(a/r), r > a,
dγ J0 (γr) = π (9.23)
0 γ 2
, r ≤ a.
Fijándonos en las condiciones (9.8) y (9.11), concluimos pues que la so-
lución buscada es la siguiente:
Z ∞
2 sen(γa)
C(r, z) = C∞ − C∞ dγ e−γz J0 (γr) . (9.24)
π 0 γ
Parra arribar a la ec. (9.1) nos preguntamos con qué velocidad alcanzan
las partı́culas al disco absorbente. Para ello, debemos calcular la tasa de
cambio de la concentración de partı́culas en la superficie del mismo, es decir,
∂C/∂z, en z = 0 y r ≤ a. Aplicando esta operación a (9.24), uno obtiene el
resultado siguiente después de calcular las integrales resultantes:
∂C(r, z) 2C∞ 1
= √ . (9.25)
∂z z=0,r≤a π a2 − r 2
El flujo superficial está dado por,
D ∂C E
Jsup = −D , (9.26)
∂z z=0
donde los paréntesis angulares significan en este caso, integrar sobre todos
los puntos (r, θ) sobre el disco de radio a, esto es,
2C∞ 2π a
 Z Z 
r
Jsup = −D dθ dr √ .
π 0 0 a2 − r2
Esta es una integral estándar y no es difı́cil llegar encontrar que su valor
es 4C∞ a, de manera que,
Jsup = −4C∞ Da. (9.27)
Por otra parte, sabemos que,
Jsup = N/τ, y C∞ = N/V.
Junto con la expresión recién obtenida, se sigue entonces que
1 4Da
= , (9.28)
τ V
donde τ es el tiempo promedio de sobrevivencia visto antes (sección 8). Su
recı́proco da la tasa de arribo k (partı́culas por unidad de tiempo) al disco
absorbente. La expresión resultante es (9.1).

33
Figura 6: Cubos conectados por un disco de radio a. El compartimento de la
derecha será etiquetado con el número 1, y el de la derecha con el número 2.

9.2. Cavidades cúbicas

Consideremos por inicio de cuentas que los contenedores son cúbicos,
véase Fig. 6. Escribiremos en seguida la ecuación de evolución del número de
partı́culas del contenedor designado por el número 1. Para esto recordaremos
que la partı́cula browniana se comporta como un caminante aleatorio (véase
la sección 2), y que por tanto al estar la partı́cula en el disco que conecta a los
volúmenes, tendrá una probabilidad de 1/2 de pasar del contenedor donde
se encontraba al volumen contiguo, y una probabilidad de 1/2 de refregar al
volumen de donde partió.
Sean N1 , N2 , y k1 , k2 , el número de partı́culas y las constantes de equilibrio
de los contenedores 1 y 2, respectivamente. La cantidad k1 N1 (t) da la fracción
de partı́culas del recipiente 1 que por unidad de tiempo alcanzan al disco,
y por analogı́a se dice lo correspondiente para el producto k2 N2 (t), como de
la primera fracción en promedio la mitad de partı́culas pasan al contenedor
2 (partı́culas que salen de V1 ), y de la segunda fracción igualmente la mitad
pasa al contenedor 1 (partı́culas que entran a V1 ), la ecuación que da la
evolución de N1 es,
dN1 1 1
= − k1 N1 (t) + k2 N2 (t). (9.29)
dt 2 2
Sustituyendo (9.1) en (9.29), esta última ecuación puede ser reescrita

34
como,
dN1 h N (t) N (t) i
2 1
= 2aD − , (9.30)
dt V2 V1
donde se supone que la constante de difusión es la misma para los dos medios.
Usando la consevación del número total de partı́culas (las paredes de los
recipientes en este caso se toman como reflejantes) N1 (t) + N2 (t) = N , uno
encuentra que,
dN1 h (N − N (t)) N (t) i
1 1
= 2aD −
dt V2 V1
hN V1 + V2 i
= 2aD − N1 (t)
V2 V1 V2
V1 + V2 h N N1 (t) i
= 2aD − . (9.31)
V2 V1 + V2 V1
Definiendo la concentración en el recipiente 1 como C1 (t) = N1 (t)/V1 y
la concentración en el equilibrio como C = N/(V1 + V2 ), la ecuación queda
como sigue,
dC1 V1 + V2
= 2aD [C − C1 (t)]. (9.32)
dt V1 V2
(0)
Considerando la condición inicial C1 (0) = C1 , y con C̃(t) = C − C1 (t),
tenemos,
dC̃(t) V1 + V2
= −2aD C̃(t), (9.33)
dt V1 V2
que tiene por solución a,
(0)
 V1 + V2 
C̃(t) = C̃ exp − 2a Dt , (9.34)
V1 V2
(0)
o bien, dado que C̃ (0) = C − C1 , obtenemos,
(0)
 V1 + V2 
C1 (t) = C + [C1 − C] exp − 2a Dt . (9.35)
V1 V2
De aquı́ notamos que la concentración sólo depende de los parámetros
geométricos del sistema y de la constante de difusión. Nótese que para tiem-
pos muy largos se recupera la concentración en el equilibrio C, tal y como es
de esperar. Se define la función de ralajación del sistema como,
C1 (t) − C  V1 + V2 
R(t) := (0) = exp − 2a Dt . (9.36)
C1 − C V1 V2

35
 La razón para hacer esta definición es que tenemos ahora una función que
da la evolución de la concentración siendo independiente del valor inicial de
la misma. También podemos definir el tiempo de relajación, dado por,
h V1 + V2 i−1
τ = − 2a D . (9.37)
V1 V2
Este parámetro da el tiempo que en promedio le toma al sistema llegar
al equilibrio.

9.3. Cavidades esféricas

Habiendo analizado la manera de proceder en este tipo de problemas,
veamos ahora la situación en que los contenedores son esféricos y las super-
fiecies de los mismos son de naturaleza absorbente (Fig. 7).Las ecuaciones
de evolución del número de partı́culas en los dos recipientes son similares a
(9.29), pero ahora debemos incluir el efecto de tener fronteras absorbentes.
Esto se hace introduciendo un tercer término de la forma −ki Ni (t), que es
el da cuenta de las partı́culas que abandonan por definitivo al sistema. Para
simplificar un tanto las cosas, tomaremos k1 = k2 = k (lo cual por (9.1), im-
plica la igualdad de los coeficientes de difusión y de los volúmenes de ambas
cavidades), con esto se tiene pues que,
dN1 1 1
= kN2 (t) − kN1 (t) − kN1 (t), (9.38)
dt 2 2

dN2 1 1
= kN1 (t) − kN2 (t) − kN2 (t), (9.39)
dt 2 2
Para resolver este sistema de ecuaciones acopladas, sumamos y restamos
a (9.38) y (9.39), para obtener,
d
(N1 (t) + N2 (t)) = −k(N1 (t) + N2 (t)), (9.40)
dt

d
(N1 (t) − N2 (t)) = 2k(N1 (t) − N2 (t)). (9.41)
dt
Resolvemos la ec. (9.40) sujeta a la condición inicial N1 (0) + N2 (0) = N0 ,
lo cual da,
N1 (t) + N2 (t) = N0 e−kt . (9.42)

36
Figura 7: Esferas conectadas por un disco de radio a. El compartimento de la
derecha será etiquetado con el número 1, y el de la derecha con el número 2.

Como las partı́culas pueden desplazarse con igual probabilidad a la iz-

quierda o derecha, en t = 0, tenemos que N1 (0) − N2 (0) = 0. Resolviendo la
ec. (9.41) con esta condición se encuentra que,

N1 (t) − N2 (t) = 0 → N1 (t) = N2 (t). (9.43)

Combinando las ecs. (9.15) y (9.16), uno encuentra que

1
Ni (t) = N1 (t) = N2 (t) = N0 e−kt , (9.44)
2
que al dividir entre el volumen de una de las esferas, da la concentración en
función del tiempo, esto es,

Ci (t) = C0 e−kt , (9.45)

con C0 = N0 /2V . Y dado que k = 4aD/V = 3aD/πR3 , con R el radio de las

esferas, tenemos la expresión siguiente,
h 3a i
Ci (t) = C0 exp − Dt . (9.46)
πR3
En este caso, el tiempo de relajación del sistema es πR3 /3aD.

37
10. Procesos de difusión en geometrı́as com-
plejas: Endocitosis. Método de los propa-
gadores
La endocitosis mediada por receptores es un mecanismo general entre las
células animales por medio del cual éstas son capaces de llevar a su inte-
rior una gran variedad de material extracelular: hormonas, glucoproteı́nas,
enzimas, toxinas, virus, etc.
El proceso se da como sigue: en un principio los ligandos llegan a sus
receptores, los cuales se encuentran esparcidos sobre la superficie de la mem-
brana celular. El complejo ası́ formado queda fijo en algunas zonas de la
membrana gracias a la acción de ciertas proteı́nas llamadas clatrinas. Pos-
teriormente, parte de la membrana se pliega sobre sı́ misma, formando una
vesı́cula con los complejos ligando–proteı́na en su interior, que pasa al seno
de la célula. Una vez en el interior de la célula, cambios de pH en el medio
provocan que los ligandos se liberen de sus receptores, a la vez que estructu-
ras tubulares se unen a la vesı́cula. Ya libres, tanto ligandos como receptores
efectúan movimiento browniano, difundiendo los receptores hacia los túbulos
anexados a la vesı́cula. En un tiempo de aproximadamente 5 a 10 minutos, en
el laboratorio se observa que el 95 % de los ligando permanecen en la vesı́cula
y que más o menos el mismo procentaje de receptores se encuentran en los
tubos, los cuales regresan a la membrana celular y liberan a los receptores,
quedando éstos en condiciones de volver a capturar nuevos ligandos.
El diámetro de las vesı́culas endocı́ticas es de aproximadamente de 200–
800 nm, mientras que los tubos que se unen a las primeras tienen un diametro
de entre 10 y 60 nm. Se sabe que del volumen total del sistema un 60 o 70 %
corresponde a la vesı́cula. Surge a partir de esta observación experimental
un punto interesante, y es que si el tiempo de relajación del sistema fuese
menor que los 5 o 10 minutos en que se reparten los ligandos y receptores de
la forma descrita en el párrafo anterior, serı́a un tanto complicado explicar
las concentraciones medidas. En efecto, cuando el sistema llega al equilibrio
la distribución de ligandos tanto en la vesı́cula como en los tubos, deberı́a ser
homogénea. De manera que si el volumen de la vesı́cula es un 60 o 70 % del
volumen total, uno esperarı́a que ese porcentaje de ligandos fuesen hallados
en el interior de las mismas, y no el 95 % que se reporta. Para evitar este
tipo de incongruencias, uno deberı́a ser capaz de mostar que el tiempo en
que relaja el sstema es mayor que 10 minutos.

38
 Para atacar tal cuestón de forma cuantitativa, podemos emplear la ecua-
ción de difusión para modelar la migración de los ligandos en el sistema
vesı́cula–tubos. Sobre esta lı́nea J. L. Lindermann y D. A. Lauffenburger8
trataron el problema pero considerando que la vesı́cula era esférica e intro-
duciendo el efecto de los tubos como un disco absorbente sobre la superficie
de la esfera. En este caso, sabemos que el tiempo de relajación viene dado
por (ver sección 9),
Vves
τ= , (10.1)
4aDves
donde Vves es el volumen de la vesı́cula, a el radio del disco por donde pue-
den escapar los ligandos, y Dves es el coeficiente de difusión de la vesı́cula.
Considerando los siguientes valores tı́picos Rves = 4 × 10−5 cm, a = 10−6 cm,
Vves = 0,7Vtotal y Dves = 10−7 cm2 /s, uno puede hacer la estimación siguiente:

0,7(2,68 × 10−13 cm3 )

τL−L = ≈ 0,47s. (10.2)
4(10−6 cm)(10−7 cm2 /s)

Si fuese correcto, la concentración de ligandos se equilibrarı́a casi de in-

mediato al iniciar el proceso, y no se podrı́a tener que despueés de 10 minutos
un 95 % de los mismos se hallase en la vesı́cula. Claramente, con las conside-
ración hechas acerca de la geometrı́a del problema lo sobresimplifican, y nos
vemos en la necesidad de considerar un modelo que no pase por alto muchos
de los detalles de la situación bajo estudio.

10.1. El método de los propagadores

Mostraremos cómo usar el método de los propagadores para calcular de el
tiempo de relajación del sistema sin tener que suponer una forma especı́fica de
la vesı́cula y con la posibilidad de introducir tantos tubos como se desee. Un
propagador, denotado por Gi (t), da la densidad de probabilidad de encontrar
una partı́cula en la cavidad i al tiempo t. Comencemos pues estableciendo
las hipótesis con las que trabajeremos. Consideraremos una vesı́cula de forma
arbitraria y volumen Vves , a la cual se le pegan tubos de radio b y longitud
L. Las constantes de difusión en la vesı́cula y en los tubos son designadas
por Dves y Dtub , respectivamente. Nuestro estudio comenzará una vez que
los ligandos se han separados de sus respectivos receptores. Supondremos que
8
Lindermann JL, Lauffenburger,Receptor/ligands sorting along the endocytic pathway,
Springer-Verlag, Berlı́n, Heidelberg (1989).

39
inicialmente todos los ligandos se hallan en el seno de la vesı́cula. Asumiremos
que se sigue una cinética de primer orden en la transición de un ligando de
la vesı́cula hacia un tubo, y por lo tanto la constante de velocidad de la
transición estará determinada por la constante de Hill (ver sección 9), es
decir,
4bDves
k= . (10.3)
Vves
En analogı́a a lo hecho en la sección anterior, podemos defirnir una fun-
ción de relajación cambiando en la foŕmula (9.9) la concentración (que co-
rresponderı́a a la vesı́cula en este caso) por el propagador Gves , y con la
concenctración en el equilibrio dada por,
Vves
Peq := lı́m Gves (t) = , (10.4)
t→∞ Vves + Vtub
quedando entonces la función de ralajación como sigue:
Gves (t) − Peq
R(t) = . (10.5)
Gves (0) − Peq
Con esta función se calcula el tiempo de relajación promedio a través de
la ecuación siguiente:
Z ∞
hτ i = dt R(t) = R̂(0), (10.6)
0

donde el acento circunflejo indica que se trata de la transformada de Laplace

de la cantidad en cuestión.9
Se infiere entonces que para develar el aspecto dinámico del sistema bajo
estudio debemos obtener una expresión para el propagador Gves (t), y es la
búsqueda en la que nos centraremos ahora. Comencemos escribiendo la ecua-
ción de evolución de Gves . Dicha ecuación, como se justificará en seguida,
está dada por,
Z t
d
Gves (t) = −kGves (t) + dσ kGves (σ)k 0 Gtub (t − σ), (10.7)
dt 0
9
La tranformada de Laplace de la función f (x), está dada por:
Z ∞
fˆ(s) = dx e−sx f (x).
0

40
con k 0 una medida de la eficiencia con la que un ligando que se encuentra en
un tubo sale del mismo hacia la vesı́cula y Gtub (t) la densidad de probabilidad
de hallar un ligando al tiempo t dentro de un tubo. Como podemos apreciar,
el cambio en la densidad de probabilidad de encontrar un ligando en el seno
de la vesı́cula al tiempo t, cambia debido a dos factores. El primero (primer
término en el lado derecho de (10.7)) da cuenta de los ligandos que alcanzan
por primera vez la entrada de uno de los tubos al tiempo t. El segundo
término tiene que ver con los ligandos que pasan de la vesı́cula a los tubos
al tiempo σ < t, difunden en los túbulos durante un intervalo de duración
t − σ, y vuelven a la vesı́cula justo al tiempo t.
El valor de k 0 es sencillo de determinar usando un argumento de flujos en
equilibrio. En efecto, con la aproximación de que el volumen de la vesácula es
inmensamente más grande que el de un tubo, y que a este último lo podemos
considerar un canal lineal, se sigue que las fracciones de ligandos en la vesı́cula
y el canal, están dadas por
Vves
fves = ≈ 1, (10.8)
Vves + Vtub

Vtub Vtub
ftub = ≈ . (10.9)
Vves + Vtub Vves
Las concentraciones respectivas se obtienen de multiplicar las fracciones
anteriores por el número total de ligandos, digamos N , y divididas por el volu-
men de la vesı́cula y la longitud del tubo (puesto que estamos aproximándolo
por un canal unidimensional). Uno obtiene lo siguiente,

N
Cves = , (10.10)
Vves

Vtub N
Ctub = . (10.11)
Vves L
Al multiplicar estas cantidades por Vves k y k 0 , respectivamente, uno en-
cuentra las corrientes de salida y entrada (referidas a si los ligandos salen o
entran de la vesı́cula), y dividiendo por la sección transversal del canal, se
llega a las expresiones para los flujos, las cuales son
kN 4bDves N
Jves = 2
= 2 , (10.12)
b b Vves

41
 k 0 N Vtub
Jtub = . (10.13)
b2 LVves
Suponiendo que los flujos se encuentran en equilibrio, podemos igualar
las dos ecuaciones anteriores y resolver para k’,

4bDves N k 0 N Vtub 0 4bDves 4Dves

= → k = = , (10.14)
b2 Vves b2 LVves Vtub πb

donde consideramos que Vtub = πb2 L.

La ecuación (10.7) ha de ser resuelta sujeta a la condición inicial mencio-
nada lı́neas arriba, a saber, que la totalidad de ligando se halle en la vesı́cula
en t = 0, esto es,
Gves (0) = 1. (10.15)
Luego, reconocemos el segundo término del miembro derecho de (10.7)
como una convolución lineal, con lo cual es natural aplicar la transformada
de Laplace y reducir tal ecuación a una algebraica. Requeriremos algunas
propiedades básica de esta transformada, a saber, la de linealidad, la tran-
formada de la derivada,
n df o
L = sfˆ(s) − f (0), (10.16)
dt
ası́ como la transformada de la convolución,

L{f ∗ g} = fˆ(s)ĝ(s). (10.17)

Aplicando pues la transformada de Lapalce a la ec. (10.7), teniendo en

cuenta las anteriores propiedades, se tiene que,

sĜves − 1 = −k Ĝves + k Ĝves k 0 Ĝtub , (10.18)

de donde se sigue que,

1
Ĝves = . (10.19)
s + k(1 − k 0 Ĝtub )

Para continuar nuestro desarrollo debemos obtener una relación explı́cita

del término k 0 Ĝtub . Para hallar dicha expresión supondremos que los ligandos
dentro de los tubos se mueven como si se tratase de una simple difusión

42
con la que ya estamos familiarizados. Reemplazaremos el tubo por un canal
unidimensional,10 Tomamos el origen del canal en el punto donde se une con
la vesı́cula. Consideramos ahora un ligando que en el tiempo t = 0 entra en
el canal. La densidad de probabilidad de encontrar a un tiempo t posterior
al ligando en el punto x, está representada por Gtub (x, t|0), densidad que
suponemos satisface la ecuación de difusión, esto es,

∂ ∂2
Gtub (x, t|0) = Dtub 2 Gtub (x, t|0). (10.20)
∂t ∂x
Analicemos qué condiciones a la frontera parecen apropiadas para este
problema. En primer lugar, la pared colocada en x = L ha de ser de natura-
leza reflejante, lo cual quiere decir que
∂
Gtub (x, t|0) = 0. (10.21)
∂x x=L

Por otro lado, no todos los ligando que alcanzan el punto x = 0 pasan a la
vesı́cula, algunos de ellos regresaran al tubo, por lo tanto supondremos que
la pared colocada allı́ tiene la caracterı́stica de ser parcialmente absorbente,
esto es,
∂
Gtub (x, t|0) = k 0 Gtub (x = 0, t|0). (10.22)
∂x x=0

Aplicamos la transformada de Laplace (sobre la variable t) a las tres

últimas expresiones, y obetenemos que,

d2
sĜtub (x, s|0) − δ(x) = Dtub Ĝtub (x, s|0), (10.23)
dx2

d
Ĝtub (x, s|0) = 0, (10.24)
dx x=L

d
Dtub Ĝtub (x, s|0) = k 0 Ĝtub (0, s|0). (10.25)
dx x=0

La función delta de Dirac δ(x), aparece en la ec. (10.23) pues tal es la

condición inicial: en t = 0 la partı́cula se halla en x = 0. La ec. (10.23) es
una ecuación diferencial ordinaria lineal de segundo orden no homogénea. Si
10
Esta una aproximación que se encuentra justificada en base a simulaciones numéricas,
donde se observa un gran acuerdo entre lo obtenido para los dos sistemas. Véase:

43
la función delta de Dirac no estuviese presente en (10.23), nos resultarı́a una
ecuación homogénea cuya solucón general es bien conocida, a saber,
r 
(h) s
Ĝtub (x, s|0) = A cosh (x − δ) , (10.26)
Dtub
y por tanto,
r r 
d (h) s s
Ĝ (x, s|0) = A senh (x − δ) , (10.27)
dx tub Dtub Dtub
con A y δ constantes a ser determinadas usando las condiciones a la frontera.
La idea básica para encontrar la solución a (10.15), será considerar un número
 > 0, y resolver la ecuación en el intervalo [, L], para después hacer  tender
a cero y recuperar la solución sobre todo el intervalo [0, L]. La solución en el
primer intervalo mencionado viene dada por (10.26), pero debemos encontrar
la condición de frontera apropiada en el nuevo extremo x = . Para ello,
introduciremos primero la noción del flujo de densidad de probabilidad, el
cual está dado por
ϕ(t) := k 0 Gtub (t) = k 0 Gtub (0, t|0), (10.28)
de manera que ϕ̂(s) = k 0 Ĝtub (0, s|0), y la condición (10.25) queda como
d
Dtub Ĝtub (x, s|0) = ϕ̂(s). (10.29)
dx x=0

Ahora bien, para obtener la condición sobre la primera derivada de Ĝtub

con respecto a x en x = , integramos la ec. (10.23) de 0 a :
Z 
dĜtub 
Dtub −s dx Ĝtub (x, s|0) = −1. (10.30)
dx 0 0

Al hacer tender  → 0 la integral en el lado izquierdo se anula. En efecto,

usando la estimación siguiente, que se prueba en los cursos elementales de
cálculo, Z
dx f (x) ≤ sup |f (x)| · longitud(γ), (10.31)
γ x∈γ

y considerando que la densidad de probabilidad está normalizada, es claro lo

anterior. Tenemos pues que,
dĜtub dĜtub
lı́m Dtub = −1 + Dtub = −1 + ϕ̂(s). (10.32)
→0 dx x= dx x=0

44
 Imponemos pues, las condiciones (10.24) y (10.32) en (10.27) para deter-
minar las constantes A y B. De considerar la condición en x = L, resulta
r r 
d s s
Ĝtub (x, s|0) =0= A senh (L − δ)
dx x=L Dtub Dtub

⇒ δ = L. (10.33)
Aplicando la condición en x = , obtenemos
r 
d p s
Dtub Ĝtub (x, s|0) = −1 + ϕ̂(s) = sDtub A senh ( − L) ,
dx x= Dtub

−1 + ϕ̂(s)
⇒A= √ hq i. (10.34)
s
sDtub senh Dtub
( − L)

Y tomando el lı́mite en que  → 0, llegamos a la solución de (10.23) que

satisface las condiciones de frontera (10.24) y (10.25):
hq i
s
1 − ϕ̂(s) cosh Dtub
(x − L)
Ĝtub (x, s|0) = √ hq i , (10.35)
sDtub senh s
L
Dtub

Empleando la definición (10.28), se sigue que

k 0 (1 − ϕ̂(s)) cosh (σ)
ϕ̂(s) = k 0 Ĝtub (0, s|0) = √ ,
sDtub senh (σ)
p
donde σ = L s/Dtub . Llegamos pues al resultado,
k0
ϕ̂(s) = √ , (10.36)
k0 + sDtub tanh (σ)
o bien, con α = k 0 L/Dtub , podemos reescribir la anterior ecuación como
σ tanh (σ)
ϕ̂(s) = 1 − . (10.37)
α + σ tanh (σ)
Sustituyendo (10.29) en (10.11), obtenemos
h σ tanh (σ) i−1
Ĝves (s) = s + k . (10.38)
α + σ tanh (σ)

45
 De esta ecuación se ve que la fracción de ligandos en la vesı́cula depende
sólo de la geometrı́a del sistema y de los coeficientes de difusión, justo como
en los ejemplos analizados en la sección anterior. Este resultado también nos
permite calcular la transformada de Laplace de la función de ralajación,
Ĝves (s) − Peq /s
R̂(s) = , (10.39)
1 − Peq
donde hemos sustituido la condición inicial Gves (t = 0) = 1 [ec. (10.15)]. No
sacaremos la transformada de Laplace inversa,únicamente nos limitaremos a
calcular el tiempo promedio de relajación, el cual viene dado por la ec. (10.6).
Para este fin, debemos llevar a cabo un desarrollo en serie de Taylor de la
función R̂(s), en torno al punto s = 0, y darnos cuenta de que
1 1
Peq = = , (10.40)
1 + kttub 1 − k ϕ̂0 (0)
donde ttub = L/k 0 = πbL/4Dves , es el tiempo promedio de vida medio de un
ligando dentro de un tubo. Primero verificaremos que la primera igualdad en
esta última ecuación es congruente con la expresión en (10.4), para lo cual
sólo debemos sustituir el valor de k [ec. (10.1)], y recordar que Vtub = πb2 L,
tenemos pues que
1 1
Peq = =   
1 + kttub 4bDves
1 + Vves πbL
4Dves
1 1 Vves
= πb2 L
= Vtub
= .
1+ Vves
1 + Vves Vves + Vtub

La segunda igualdad en (10.40) se comprueba haciendo un desarrollo en

serie de Taylor en torno al punto s = 0 de la derivada ϕ̂0 (s). Tal desarrollo
lo mostraremos más adelante, por ahora daremos por cierta la ec. (10.40)
y proseguiremos en el cálculo del tiempo promedio de relajación. Hacien-
do un desarrollo en serie de Taylor de ϕ̂(s), alrededor del punto s = 0, y
quedándonos con los primeros tres términos, obtenemos que,
1
ϕ̂(s) ≈ ϕ̂(0) + ϕ̂0 (0)s + ϕ̂00 (0)s2 . (10.41)
2
Sustiuyendo (10.41) en (10.19), se obtinene que,
1
Ĝves (s) ≈ . (10.42)
s− ks(ϕ̂0 (0) + 21 ϕ̂00 (0)s)

46
 E introduciendo (10.42) junto con (10.40) en (10.39), se obtiene lo si-
guiente,
1
s−ks(ϕ̂0 (0)+ 12 ϕ̂00 (0)s)
− s(1−k1ϕ̂0 (0))
R̂(s) =
1 − 1−k1ϕ̂0 (0)
1−kϕ̂0 (0)
s−ks(ϕ̂0 (0)+ 12 ϕ̂00 (0)s)
− 1s
=
−k ϕ̂0 (0)
s − k ϕ̂0 (0)s − s + sk(ϕ̂0 (0) + 21 ϕ̂00 (0)s)
=
−k ϕ̂0 (0)s2 [1 − k(ϕ̂0 (0) + 21 ϕ̂00 (0)s)]
1 00
2
ϕ̂ (0)
= ,
−ϕ̂ (0)[1 − k ϕ̂0 (0) − k2 ϕ̂00 (0)s]
0

y con s = 0, obtenemos,
1 00
2
ϕ̂ (0)
R̂(0) = . (10.43)
−ϕ̂0 (0)[1 − k ϕ̂0 (0)]
Para conocer los valores de las derivadas que aparecen en la ecuación an-
terior, hacemos un desarrollo en serie de Taylor de (10.37) y de sus derivadas,
en torno al punto s = 0. Dichos desarrollorros se muestran a continuación:

L3 L2
 
L
+ 02 s2 + O s3 ,


ϕ̂(s) = 1 − 0 s + 0
(10.44)
k 3Dtub k k

L 2L2 (3Dtub + k 0 L)
ϕ̂0 (s) = − s + O s3/2 ,


0
+ 02
(10.45)
k 3Dtub k

2L2 (3Dtub + k 0 L) 2 [L3 (15Dtub

2
+ 10Dtub k 0 L + 2k 02 L2 )]
ϕ̂00 (s) = 3/2


− 2
s+O s .
3Dtub k 02 5 (Dtub k 03 )
(10.46)
De aquı́ se sigue que,
L
ϕ̂0 (0) = − = −ttub , (10.47)
k0
y,
2L2 (3Dtub + k 0 L) L2
 
00
ϕ̂ (0) = = 2ttub ttub + . (10.48)
3Dtub k 02 3Dtub

47
 Al sustituir (10.47) y (10.48) en (10.43), encontramos que
L2
ttub + 3Dtub
hτ i = , (10.49)
1 + kttub
o bien,
L2
 
Vves πbL
hτ i = + . (10.50)
Vves + Vtub 4Dves 3Dtub
Nuevamente resalta la dependencia únicamente en los parámetros geométri-
cos del sistema y en las constantes de difusión. Para comparar el resultado
obtenido con el método de los propagadores cuando se toma en cuenta la
presencia de los tubos con respecto a no considerarlos, multiplicamos por k
la ecuación anterior y usamos los mismoas valores numéricos que antes para
obtener,
khτ i ≈ 3260. (10.51)
Es decir, que el tiempo promedio de relajación calculado aquı́ es unas 3260
veces mayor que el obtenido antes [ec. (10.2)]. La principal consecuencia de
esto, es que después de diez minutos el sistema está muy lejos del equilibrio
y por ende los ligandos no han alcanzado a difundir a los tubos, razón por la
cual la gran mayorı́a de ellos permanece en la vesı́cula.

48
 Difusión en sistemas confinados simétricos en
2D y 3D.
Reporte 2.

Roberto Verdel Aranda

3 de junio de 2014

1. Preliminares
En la primera parte del presente proyecto, introdujimos el estudio de la
difusión bajo confinamiento con geometrı́as simples, e inclusive tratamos un
caso donde la geometrı́a hace del problema (endocitosis) uno no trivial. En
aquella ocasión, empleamos el método de los propagadores, con el que fuimos
capaces de obtener el tiempo promedio del primer arribo a un cierto dominio,
lo que permitió esclarecer la parte dinámica del problema. Se estudió tam-
bién el problema de difusión entre dos cavidades interconectadas entre sı́. La
estrategia seguida entonces fue la de buscar cómo cambia la concentración
en el tiempo, en lugar de tratar de determinar a ésta como una función del
punto y del tiempo.
En esta segunda parte del proyecto, continuaremos nuestro estudio de
la difusión en sistemas confinados, centrándonos ahora en sistemas en dos
y tres dimensiones que presenten simetrı́as fáciles de identificar. La traza
que seguiremos en esta ocasión será distinta a las comentadas en el párrafo
anterior. La idea aquı́ será introducir el potencial inducido por el confina-
miento en términos de un coeficiente de difusión efectivo en la ecuación de
Smoluchowski.
Una de las principales motivaciones para llevar a cabo estudios de esta
naturaleza, es la inmensa cantidad de aplicaciones que pueden encontrarse de
la difusión en sistemas confinados en múltiples contextos dentro de la fı́sica,
la quı́mica, la biologı́a y en la nanotecnologı́a. Por razones de espacio, no

1
desglosaremos la enorme lista de referencias con las aplicaciones a las que
hemos hecho alusión. En su lugar, referimos al lector a consultar la Ref.[?],
donde se presenta un listado detallado de tales aplicaciones.

1.1. Ecuación de Fick–Jacobs–Zwanzig y coeficientes

de difusión efectivos
En el primer reporte de este trabajo, vimos que el transporte por difusión
pasiva está descrito por la segunda ley de Fick o ecuación de difusión:

∂C(~r, t)
= D0 ∇2 C(~r, t), (1.1)
∂t
donde D0 es la constante de difusión y C(~r, t) es la concentración de partı́culas
como función del punto y del tiempo.
En su momento se comentó la dificultad que acompaña la solución de
(1.1) con condiciones de frontera arbitrarias al igual que la condición inicial.
La principal lı́nea seguida dentro de los esfuerzos teóricos para superar esta
complicación, ha sido la de hacer una descripción simplificada mediante una
reducción dimensional del problema en cuestión. Esto es, se trata de buscar
un problema equivalente al original en donde la ecuación de evolución esté en
términos de únicamente cantidades unidimensionales. En varios problemas
cuasi-unidimensionales, es la dirección longitudinal del sistema la que resulta
de mayor interés. Las coordenadas correspondientes a la parte transversal no
se toman en cuenta bajo el supuesto de que las partı́culas llegan al equilibrio
casi instantáneamente en dicha dirección.
La primera contribución en esta manera de proceder se atribuye a Mer-
kel H. Jacobs, quien en 1935 hace una publicación donde aparece una de-
ducción heurı́stica de la ecuación actualmente conocida como ecuación de
Fick–Jacobs:  
∂C(x, t) ∂ ∂ C(x, t)
= D0 w(x) , (1.2)
∂t ∂x ∂x w(x)
donde en este caso, C(x, t) representa la concentración lineal efectiva en la
dirección del eje coordenado x, y w(x) es el perfil de la sección transversal
del sistema donde ocurre la difusión. La deducción rigurosa de esta ecuación
fue hecha por Robert Zwanzig en 1992. Para este fin, Zwanzig se valió de la
ecuación de Smoluchowsky:
 
∂C(x, t) ∂ −βU (x) ∂ βU (x)
= D0 e e C(x, t) , (1.3)
∂t ∂x ∂x

2
con β = 1/kB T , donde kB es la constante de Boltzmann y T la temperatura
absoluta del sistema. Esta ecuación describe la difusión de partı́culas en una
dimensión bajo la influencia del potencial externo U (x).
De comparar las ecuaciones (1.2) y (1.3), uno puede definir el potencial
entrópico:
1
U (x) := − ln[w(x)]. (1.4)
β
Bajo esta perspectiva, vemos que la descripción efectiva del problema de
confinamiento se lleva a cabo considerando un problema “libre”(en el sentido
de la limitación espacial), con el precio de tener que añadir un potencial que
lleva consigo la información de la geometrı́a del sistema, i.e. del confinamien-
to, a través de la función w(x).
Zwanzig se percató además de que la ecuación de Fick–Jacobs no era
la más adecuada para estudiar difusión en muchos casos, de manera que la
modificó introduciendo la noción de coeficientes de difusión efectivos, D(x),
los cuales dependen de la coordenada espacial x. La ecuación a la que Zwanzig
llegó, conocida hoy en dı́a bajo el tı́tulo de ecuación de Fick–Jacobs–Zwanzig,
es la siguiente,
 
∂ ∂ ∂ C(x, t)
C(x, t) = D0 D(x)w(x) . (1.5)
∂t ∂x ∂x w(x)
Después de estos trabajos, muchos fı́sicos y matemáticos se han dado a la
tarea de encontrar una forma funcional de los coeficientes de difusión efecti-
vos. En 2005, Pavol Kalinay y Jerome K. Percus, establecieron por primera
vez un método riguroso por medio del cual fue posible obtener expresiones
de estos coeficientes para sistemas confinados en geometrı́as con ciertos ele-
mentos de simetrı́a. Los detalles de este método, se discuten en la siguiente
sección.

2. Método de Kalinay y Percus

Como se mencionó con anterioridad, describiremos aquı́ un método rigu-
roso, desde el punto de vista de las matemáticas; a través del cual, proyec-
taremos la ecuación de difusión sobre la dirección longitudinal de un canal
simétrico estrecho y de sección transversal variable. Esto lo haremos tanto en
dos como en tres dimensiones. El método que a continuación presentaremos
es un mapeo basado en la teorı́a de las perturbaciones, y fue introducido
originalmente por Pavol Kalinay y Jerome K. Percus.

3
 y

J2 d Hx, wHxL, tL
w(x)

J2 d Hx, 0, tL
x
0 L

Figura 1: Canal plano de longitud L con una pared coincidente con el eje
x, y la otra de forma arbitraria w(x). El flujo es paralelo a las paredes del
canal, las cuales son de naturaleza totalmente reflejante.

2.1. Canal bidimensional con una pared plana y otra

de forma arbitraria
Comenzaremos estudiando un canal plano, i. e. bidimensional. Este canal
es de longitud L y está constituido por dos paredes, una a lo largo del eje x y
la segunda de forma arbitraria w(x). La representación esquemática de este
canal se muestra en la Fig. 1. Las dos paredes son de naturaleza reflejante,
esto es, el flujo a través de ellas es nulo. Supondremos que el flujo en el
interior del canal es paralelo a las paredes del mismo. Por otra parte, en
los bordes de éste, x = 0, L, podemos imponer condiciones a la frontera de
cualquier tipo (Dirichlet, Neumann o Robin).
La densidad bidimensional de partı́culas σ(x, y, t) (número de partı́culas
por unidad de área), satisface la ecuación de difusión, esto es,

∂σ ∂ 2σ ∂ 2σ
= Dx 2 + Dy 2 , (2.1)
∂t ∂x ∂y
donde Dx y Dy son los coeficientes de difusión en las direcciones longitudinal y

4
transversal, respectivamente. Cuando se trate de un medio isótropo, entonces
Dx = Dy .
Definimos la proyección unidimensional de la densidad, que llamaremos
concentración reducida, como
Z w(x)
P2d (x, t) := dy σ(x, y, t). (2.2)
0

Integramos ahora la ec. (2.1) sobre la dirección transversal y en todo el

ancho del canal, o sea desde y = 0 hasta y = w(x):
w(x) w(x) w(x)
∂ 2σ ∂ 2σ
Z Z Z
∂σ
dy = dy Dx 2 + dy Dy . (2.3)
0 ∂t 0 ∂x 0 ∂y 2
La integral en el miembro izquierdo de la expresión anterior se lleva a cabo
de forma inmediata. En tanto que para la primera integral en el miembro
derecho, hacemos uso de la regla de Leibnitz para derivar integrales:
Z ψ(x) Z ψ(x)
d ∂F (x, y) dψ(x) dφ(x)
dy F (x, y) = + dy F (x, ψ(x))− F (x, φ(x)).
dx φ(x) φ(x) ∂x dx dx
(2.4)
Usando este resultado tenemos que
Z w(x)  
∂ ∂σ(x, y, t)
dy Dx =
0 ∂x ∂x
( Z )
w(x)
∂ ∂σ(x, y, t) ∂
= Dx dy − w0 (x) σ(x, w(x), t)
∂x 0 ∂x ∂x
( Z w(x) )
∂2 ∂ 0 ∂
= Dx dy σ(x, y, t) − [w (x)σ(x, w(x), t)] − w0 (x) σ(x, w(x), t) ,
∂x2 0 ∂x ∂x

donde w0 (x) =dw/dx.

La segunda integral en el miembro derecho de (2.4), se integra directa-
mente usando el teorema fundamental del cálculo, dando lo siguiente,
w(x)
∂2
Z
∂ y=w(x)
dy Dy 2
σ(x, y, t) = Dy σ(x, y, t) .
0 ∂y ∂y y=0

5
 De manera que uno obtiene lo siguiente,
(
∂ ∂2 ∂ 0
P2d (x, t) = Dx P 2d (x, t) − [w (x)σ(x, w(x), t)] +
∂t ∂x2 ∂x
)
∂ ∂ y=w(x)
− w0 (x) σ(x, w(x), t) + Dy σ(x, y, t) ,
∂x ∂y y=0

(2.5)

donde hemos hecho uso de la definición (2.2).

Ahora bien, está ecuación está sujeta a las condiciones de frontera sobre
las paredes del canal, i. e. que el flujo sea paralelo al vector unitario en la di-
rección de éstas. Determinemos la expresión matemática de tales condiciones.
Según la primera ley de Fick, se tiene que,
∂σ(x, y, t) ∂σ(x, y, t)
J~2d (x, y, t) = −Dx x̂ − Dy ŷ, (2.6)
∂x ∂y
donde x̂ y ŷ son los vectores unitarios en la dirección de los ejes coordenados
x y y, respectivamente. Por como estamos ubicando las paredes (ver Fig.
1), un vector unitario en la dirección de la pared horizontal es v̂inf = x̂,
mientras que un vector unitario en la paralelo a la pared superior es v̂sup =
(x̂ + w0 (x)ŷ)/(1 + [w0 (x)]2 )1/2 , que se sigue de la interpretación geométrica de
la derivada como la pendiente de la tangente a la curva en cuestión. Entonces,
las condiciones que deseamos imponer adoptan la forma siguiente:
1 0
v̂inf × J~2d = ẑ = 0,
−Dx ∂x −Dy ∂σ
∂σ
∂y

y,
1 1 w0 (x)
v̂sup × J~2d = p ẑ = 0,
1 + [w0 (x)]2 −Dx ∂x −Dy ∂σ
∂σ
∂y

o bien, al realizar las operaciones indicadas por los determinantes, llegamos

a las siguientes expresiones:
∂
−Dy σ(x, y, t) = 0, (2.7)
∂y y=0

y,
∂ ∂
Dy σ(x, y, t) = w0 (x)Dx σ(x, y, t) . (2.8)
∂y y=w(x) ∂x y=w(x)

6
 Tomando en cuenta estas condiciones en (2.5), tenemos,
(
∂ ∂2 ∂ 0
P2d (x, t) = Dx 2
P2d (x, t) − [w (x)σ(x, w(x), t)] +
∂t ∂x ∂x
)
∂ ∂
− w0 (x) σ(x, w(x), t) − Dy σ(x, 0, t) +
∂x ∂y
| {z }
0
∂
+ Dy σ(x, w(x), t),
∂y
| {z }
∂
w0 (x)Dx ∂x σ(x,w(x),t)

es decir,
∂2
 
∂ ∂ 0
P2d (x, t) = Dx P2d (x, t) − [w (x)σ(x, w(x), t)] . (2.9)
∂t ∂x2 ∂x
Si ahora uno considera una tasa de difusión infinita en la dirección trans-
versal y, esto es, que en dicha dirección la concentración equilibra muy rápido,
que es la definición de un canal estrecho, entonces de (2.2) se sigue que
P2d (x, t)
σ(x, y, t) = . (2.10)
w(x)
Introduciendo (2.10) en (2.9), llegamos a
w0 (x)
 2  
∂ ∂ ∂
P2d (x, t) = Dx P2d (x, t) − P2d (x, t) ,
∂t ∂x2 ∂x w(x)
o bien,
w0 (x)
 
∂ ∂ ∂P2d (x, t)
P2d (x, t) = Dx − P2d (x, t) ,
∂t ∂x ∂x w(x)
donde el término entre los paréntesis cuadrados puede ser reescrito de la
siguiente manera
∂P2d (x, t) w0 (x) ∂ P2d (x, t)
− P2d (x, t) = w(x) ,
∂x w(x) ∂x w(x)
como puede verificarse fácilmente. La expresión con la que terminamos es
pues,  
∂ ∂ ∂ P2d (x, t)
P2d (x, t) = Dx w(x) , (2.11)
∂t ∂x ∂x w(x)

7
que al poner Dx = D0 , se convierte en la ecuación de Fick–Jacobs [ec. (1.2)].
Describiremos a continuación la parte central del esquema propuesto por
Kalinay y Percus. En tal esquema la ec. (2.11) representa la correción a orden
cero de la ecuación de difusión unidimensional proyectada, donde el paráme-
tro del desarrollo en series es el cociente entre los coeficientes de difusión
longitudinal y transversal λ = Dx /Dy . Enseguida mostraremos un enfoque
de recurrencia que permita obtener correcciones a orden superior de la ecua-
ción de Fick–Jacobs. Para esto, escribimos la densidad bidimensional como
una serie tal y como lo dicta la teorı́a de perturbaciones, esto es,
∞
X
σ(x, y, t) = λj σj (x, y, t), (2.12)
j=0

donde σj (x, y, t) tiene la forma de algún operador actuando sobre P2d (x, t)/w(x),
de modo tal que (2.9) sea consistente.