CCEl Club de la Ciencia
Neurociencias · Laboratorios

Neurociencias · Unidad 02

Laboratorio del físico

Del movimiento aleatorio a la ecuación de difusión

La membrana neuronal: una frontera dinámica

El laboratorio del físico

Del movimiento aleatorio a la ecuación de difusión

Este laboratorio es una lectura opcional. Puedes estudiar sus primeras secciones y regresar más adelante a las deducciones y los métodos numéricos más avanzados; no es necesario dominar todo el formalismo matemático para continuar con el texto principal de la unidad.

Nota sobre la consulta sin conexión: las ecuaciones de este documento se muestran mediante MathJax, cargado desde Internet. Para visualizarlas correctamente, abre el HTML con conexión. Los archivos de Python y la animación que aparecen en la sección de descargas pueden guardarse para utilizarlos por separado sin conexión, una vez instaladas las dependencias correspondientes.

Pregunta científica: ¿Cómo puede el movimiento aleatorio de un número muy grande de partículas producir un flujo neto de materia y cómo podemos predecir la evolución de su concentración?

En el texto principal hemos estudiado la difusión como uno de los mecanismos fundamentales del transporte de sustancias. Sabemos que las partículas se encuentran en movimiento continuo y que, cuando existen diferencias de concentración, puede producirse un transporte neto de materia.

Este hecho nos conduce a una pregunta central. Aunque la trayectoria de cada partícula individual sea irregular y no pueda predecirse en detalle, el comportamiento colectivo de un gran número de partículas sí puede describirse matemáticamente en términos estadísticos. Nuestro objetivo será entender cómo, a partir de movimientos microscópicos desordenados, emergen regularidades macroscópicas que pueden expresarse mediante leyes físicas.

Para responder, construiremos un modelo que relaciona el movimiento microscópico de las partículas con la evolución macroscópica de la concentración. Introduciremos el concepto de flujo, estudiaremos la primera ley de Fick y utilizaremos el principio de conservación de la materia para obtener la ecuación de difusión.

También examinaremos la relación entre este modelo continuo y una descripción probabilística basada en caminatas aleatorias. Esta conexión nos permitirá comprender cómo una ecuación diferencial determinista puede surgir como descripción efectiva del comportamiento estadístico de un gran número de partículas.

Finalmente, desarrollaremos una simulación numérica en Python para estudiar la evolución temporal de una distribución de concentraciones, comprobar la conservación de la materia y examinar las limitaciones del modelo.

1. El problema físico

Imaginemos un recipiente alargado que contiene agua. Inicialmente, una sustancia disuelta se encuentra concentrada en una pequeña región del recipiente, mientras que su concentración es menor en las regiones restantes.

Con el paso del tiempo, las partículas de la sustancia se dispersan por el líquido y la distribución de concentraciones tiende a hacerse más uniforme. Este proceso ocurre sin que sea necesario imponer una corriente macroscópica al líquido: es consecuencia del movimiento térmico de las partículas y recibe el nombre de difusión.

Aunque cada partícula experimenta un movimiento irregular, existe una tendencia estadística al transporte neto de materia desde las regiones de mayor concentración hacia las de menor concentración. Esta tendencia no implica que todas las partículas se desplacen en la misma dirección. En realidad, las partículas atraviesan continuamente cualquier superficie imaginaria en ambos sentidos, pero, cuando existe un gradiente de concentración, los flujos en sentidos opuestos no se compensan.

Para describir cuantitativamente este fenómeno necesitamos responder dos preguntas: ¿cuánta materia atraviesa una superficie por unidad de tiempo y cómo depende esa cantidad de las diferencias de concentración?

Estas preguntas nos conducirán al concepto de flujo de materia y a la primera ley de Fick.

2. Concentración y flujo de materia

Para estudiar cuantitativamente la difusión necesitamos construir un modelo que describa cómo se distribuye una sustancia en el espacio y cómo cambia esa distribución con el tiempo.

Comenzaremos por establecer algunas hipótesis que simplifican el problema:

Estas hipótesis no pretenden describir toda la complejidad de un sistema biológico. Su finalidad es aislar el fenómeno que deseamos estudiar y permitirnos identificar las relaciones físicas que lo gobiernan.

2.1. La concentración

La primera magnitud que necesitamos es la concentración, que nos indica cuánta sustancia se encuentra en una determinada región del espacio.

Representaremos la concentración mediante la función

\[c=c(x,t),\]

donde \(x\) es la posición, expresada en metros, y \(t\) es el tiempo, expresado en segundos.

Utilizaremos la concentración molar, entendida como la cantidad de sustancia por unidad de volumen. Sus unidades en el Sistema Internacional son

\[[c]=\mathrm{mol/m^3}.\]

La función \(c(x,t)\) contiene información tanto espacial como temporal. Nos permite conocer cómo se distribuye la sustancia en un instante determinado y cómo evoluciona esa distribución a medida que transcurre el tiempo.

Por ejemplo, si representamos gráficamente \(c(x,t)\) para un instante fijo, obtenemos un perfil de concentración. Su forma nos permite identificar las regiones donde la sustancia está más concentrada y aquellas donde su concentración es menor.

2.2. El flujo de materia

Conocer la concentración no es suficiente para describir el transporte. También necesitamos una magnitud que nos indique cuánta sustancia atraviesa una superficie durante un intervalo de tiempo.

Para ello introducimos la densidad de flujo molar, que representaremos mediante \(J(x,t)\).

En una dimensión, el flujo se define con respecto a una superficie perpendicular al eje \(x\). Si una cantidad neta de sustancia \(\Delta n\) atraviesa una superficie de área \(A\) durante un intervalo de tiempo \(\Delta t\), el flujo medio correspondiente es

\[J=\frac{\Delta n}{A\,\Delta t}.\]

En esta expresión, \(\Delta n\) representa la cantidad neta de sustancia que atraviesa la superficie en el sentido elegido como positivo, expresada en moles.

Por consiguiente, las unidades del flujo molar son

\[[J]=\mathrm{\frac{mol}{m^2\,s}}.\]

Cuando el transporte varía con la posición o con el tiempo, consideramos intervalos suficientemente pequeños para definir el flujo local \(J(x,t)\).

Es importante distinguir entre el movimiento de las partículas individuales y el flujo neto de materia. Aunque las partículas atraviesen continuamente una superficie en ambas direcciones, el flujo neto depende de la diferencia entre las cantidades transportadas en uno y otro sentido.

Adoptaremos el convenio de que \(J\) es positivo cuando el transporte neto ocurre en el sentido positivo del eje \(x\), y negativo cuando ocurre en el sentido contrario.

2.3. ¿Cómo se relacionan la concentración y el flujo?

Hemos introducido dos magnitudes que describen aspectos diferentes del mismo fenómeno. La concentración caracteriza la distribución espacial de la sustancia, mientras que el flujo describe su transporte.

Podemos anticipar que existe una relación entre ellas: si la concentración es diferente en dos regiones, puede producirse un transporte neto de materia entre ambas.

Sin embargo, para construir un modelo cuantitativo necesitamos establecer cómo depende el flujo de las variaciones espaciales de la concentración.

Esta relación constituye el contenido de la primera ley de Fick, que estudiaremos a continuación.

3. La primera ley de Fick

En la sección anterior introdujimos dos magnitudes fundamentales: la concentración \(c(x,t)\), que describe cómo se distribuye una sustancia, y el flujo \(J(x,t)\), que cuantifica su transporte.

Ahora necesitamos establecer una relación matemática entre ambas.

Sabemos que, cuando una sustancia se encuentra distribuida de manera desigual en un medio, existe una tendencia al transporte neto de materia desde las regiones de mayor concentración hacia las de menor concentración. La primera ley de Fick permite describir cuantitativamente este comportamiento.

3.1. El gradiente de concentración

Consideremos dos regiones muy próximas de un recipiente, situadas en las posiciones \(x\) y \(x+\Delta x\). Si sus concentraciones son diferentes, podemos calcular la variación de concentración por unidad de longitud:

\[\frac{\Delta c}{\Delta x}=\frac{c(x+\Delta x,t)-c(x,t)}{\Delta x}.\]

Esta expresión nos indica cuánto cambia la concentración cuando nos desplazamos una cierta distancia en el medio.

Para describir la variación local de la concentración, consideramos el límite en el que la separación entre ambas posiciones tiende a cero:

\[\frac{\partial c}{\partial x}=\lim_{\Delta x\to 0}\frac{c(x+\Delta x,t)-c(x,t)}{\Delta x}.\]

Esta derivada espacial es la componente del gradiente de concentración en la dirección \(x\).

Su signo proporciona información sobre la distribución espacial de la sustancia. Si es positivo, la concentración aumenta al avanzar en el sentido positivo del eje \(x\); si es negativo, disminuye. Su magnitud indica qué tan rápidamente cambia la concentración con la posición.

En el Sistema Internacional, sus unidades son

\[\left[\frac{\partial c}{\partial x}\right]=\mathrm{\frac{mol}{m^4}}.\]

El gradiente es, por tanto, una magnitud que caracteriza la variación espacial de la concentración, no una medida directa de la cantidad de materia transportada.

3.2. La relación entre el flujo y el gradiente

Para numerosos sistemas, cuando el medio es homogéneo y las condiciones permanecen suficientemente próximas al equilibrio local, el flujo difusivo es proporcional al gradiente de concentración.

Esta relación constituye la primera ley de Fick:

\[\boxed{J(x,t)=-D\frac{\partial c(x,t)}{\partial x}}\]

En esta ecuación, \(J(x,t)\) es la densidad de flujo molar, \(c(x,t)\) es la concentración y \(D\) es el coeficiente de difusión.

El coeficiente \(D\) es una propiedad del sistema formado por la sustancia que se difunde y el medio en el que se encuentra. Su valor depende, entre otros factores, de la temperatura, de las propiedades del medio y de las características de las partículas.

Para un medio homogéneo e isótropo, \(D\) es una magnitud escalar positiva. Sus unidades se obtienen directamente de la primera ley de Fick:

\[[D]=\mathrm{\frac{m^2}{s}}.\]

Cuanto mayor sea el coeficiente de difusión, mayor será la magnitud del flujo para un mismo gradiente de concentración, siempre que se mantengan las condiciones en las que la ley es aplicable.

3.3. ¿Por qué aparece un signo negativo?

El signo negativo de la primera ley de Fick tiene un significado físico fundamental: indica que el flujo difusivo neto se dirige en sentido contrario al gradiente de concentración.

Para comprenderlo, imaginemos que la concentración disminuye conforme avanzamos en el sentido positivo del eje \(x\). En ese caso,

\[\frac{\partial c}{\partial x}<0.\]

Como el coeficiente de difusión es positivo, la primera ley de Fick establece que

\[J>0.\]

Esto significa que el transporte neto ocurre en el sentido positivo del eje \(x\), desde la región de mayor concentración hacia la de menor concentración.

Por el contrario, si la concentración aumenta con \(x\), el gradiente es positivo y el flujo es negativo. El transporte neto ocurre entonces en el sentido contrario al eje \(x\).

cxMayor concentraciónMenorFlujo JGradiente ∇cJ = −D ∂c/∂xD > 0
Figura 1. Representación esquemática de un perfil de concentración decreciente. El gradiente apunta hacia donde aumenta la concentración, mientras que el flujo difusivo neto tiene el sentido contrario. Se supone un coeficiente de difusión positivo y constante.

Si la concentración es uniforme, el gradiente desaparece y, de acuerdo con la primera ley de Fick, el flujo difusivo neto es nulo. Esto no significa que las partículas dejen de moverse, sino que los transportes en sentidos opuestos se compensan estadísticamente.

El signo negativo, por tanto, no es una simple convención algebraica. Expresa una propiedad fundamental de la difusión ordinaria: en las condiciones descritas por este modelo, el transporte neto tiende a reducir las diferencias espaciales de concentración.

3.4. Alcance y limitaciones de la primera ley de Fick

La primera ley de Fick es una relación constitutiva: establece cómo responde el flujo a un gradiente de concentración bajo determinadas condiciones físicas.

No es una ley universal que describa cualquier forma de transporte de materia. En sistemas donde existen interacciones importantes entre partículas, medios heterogéneos, transporte activo o fuerzas externas, pueden ser necesarias relaciones más generales.

Esta distinción resulta especialmente importante en los sistemas biológicos. Los iones que atraviesan una membrana neuronal, por ejemplo, no responden únicamente a diferencias de concentración, sino también a las fuerzas eléctricas que actúan sobre ellos.

El modelo que hemos construido describe exclusivamente el transporte difusivo. Más adelante ampliaremos esta descripción para incorporar los efectos eléctricos y comprender el transporte de partículas cargadas a través de las membranas celulares.

Por ahora, la primera ley de Fick nos permite calcular el flujo a partir de la distribución de concentraciones. Sin embargo, todavía no nos indica cómo cambia esa distribución con el tiempo.

Para responder a esta última pregunta necesitaremos introducir un principio físico adicional: la conservación de la materia.

4. Conservación de la materia y ecuación de difusión

La primera ley de Fick nos permite calcular el flujo de una sustancia cuando conocemos su distribución de concentraciones. Sin embargo, nuestro objetivo es más ambicioso: queremos determinar cómo cambia esa distribución a medida que transcurre el tiempo.

Para conseguirlo necesitamos incorporar un principio físico fundamental: la conservación de la materia.

Si una región del espacio no contiene fuentes ni sumideros, la cantidad de una sustancia que se encuentra en su interior solo puede cambiar mediante el transporte a través de sus fronteras. La sustancia no aparece ni desaparece espontáneamente.

Este principio nos permitirá obtener una ecuación que relaciona la variación temporal de la concentración con el flujo de materia.

4.1. Balance de materia en un elemento de volumen

Consideremos nuevamente nuestro recipiente alargado y seleccionemos una pequeña región comprendida entre las posiciones \(x\) y \(x+\Delta x\).

Supondremos que el recipiente tiene una sección transversal constante de área \(A\). El volumen de la región considerada será, por tanto,

\[\Delta V=A\,\Delta x.\]

La cantidad de sustancia contenida en ese volumen, expresada en moles, es aproximadamente

\[n(x,t)\simeq c(x,t)\,A\,\Delta x,\]

cuando el intervalo espacial es suficientemente pequeño para considerar que la concentración es prácticamente uniforme en su interior.

AJ(x,t)AJ(x+Δx,t)Elemento de volumenΔV = A Δxxx + Δxx
Figura 2. Esquema de un elemento de volumen de sección transversal \(A\) y longitud \(\Delta x\). El flujo molar entra por la superficie izquierda y sale por la derecha. Las flechas representan los flujos netos considerados positivos en la dirección del eje \(x\); no representan las trayectorias individuales de las partículas.

El flujo molar que atraviesa la superficie izquierda del elemento es \(J(x,t)\), mientras que el correspondiente a la superficie derecha es \(J(x+\Delta x,t)\).

Recordemos que el flujo se expresa como cantidad de sustancia por unidad de área y por unidad de tiempo. Por tanto, la cantidad de sustancia que entra por unidad de tiempo es

\[\text{Entrada}=A\,J(x,t),\]

y la que sale es

\[\text{Salida}=A\,J(x+\Delta x,t).\]

En estas expresiones utilizamos flujos con signo, de acuerdo con el convenio establecido en la sección anterior. Si alguno de ellos es negativo, representa un transporte en sentido contrario al eje \(x\).

Como no existen fuentes ni sumideros, el principio de conservación establece que

\[\boxed{\text{Acumulación}=\text{Entrada}-\text{Salida}}\]

Este balance es el punto de partida para obtener nuestra ecuación de evolución.

4.2. La ecuación de continuidad

La variación temporal de la cantidad de sustancia contenida en nuestro elemento de volumen debe ser igual a la diferencia entre lo que entra y lo que sale:

\[\frac{\partial n}{\partial t}=AJ(x,t)-AJ(x+\Delta x,t).\]

Sustituyendo la expresión aproximada de \(n\), obtenemos

\[A\Delta x\,\frac{\partial c}{\partial t}=A\left[J(x,t)-J(x+\Delta x,t)\right].\]

Dividimos ambos miembros entre el volumen del elemento, \(A\Delta x\):

\[\frac{\partial c}{\partial t}=-\frac{J(x+\Delta x,t)-J(x,t)}{\Delta x}.\]

Finalmente, al tomar el límite \(\Delta x\to 0\), encontramos

\[\boxed{\frac{\partial c}{\partial t}=-\frac{\partial J}{\partial x}}\]

Esta expresión recibe el nombre de ecuación de continuidad unidimensional.

Su significado físico es sencillo: la concentración en una región aumenta cuando entra más materia de la que sale y disminuye cuando sale más materia de la que entra.

Si el flujo es uniforme en el espacio, su derivada espacial es nula. En consecuencia, la concentración local no cambia con el tiempo, aunque pueda existir transporte de materia a través de la región.

Es importante observar que la ecuación de continuidad no depende de que el transporte ocurra por difusión. Se trata de una expresión local del principio de conservación de la materia, válida también para otros mecanismos de transporte, siempre que se utilice el flujo total correspondiente y no existan fuentes ni sumideros.

4.3. Obtención de la ecuación de difusión

Disponemos ahora de dos relaciones matemáticas que describen aspectos diferentes de nuestro sistema.

Por una parte, la primera ley de Fick relaciona el flujo difusivo con el gradiente de concentración:

\[J=-D\frac{\partial c}{\partial x}.\]

Por otra, la conservación de la materia exige que

\[\frac{\partial c}{\partial t}=-\frac{\partial J}{\partial x}.\]

Podemos combinar ambas expresiones sustituyendo la primera en la segunda:

\[\frac{\partial c}{\partial t}=-\frac{\partial}{\partial x}\left(-D\frac{\partial c}{\partial x}\right).\]

Si el coeficiente de difusión es constante, podemos extraerlo de la derivada espacial y escribir

\[\boxed{\frac{\partial c}{\partial t}=D\frac{\partial^2c}{\partial x^2}}\]

Esta es la ecuación de difusión unidimensional, también conocida como la segunda ley de Fick.

Hemos obtenido una ecuación diferencial parcial que relaciona la evolución temporal de la concentración con su segunda derivada espacial.

La primera derivada espacial representa la variación local de la concentración, mientras que la segunda derivada describe cómo cambia esa pendiente de un punto a otro. En términos geométricos, está relacionada con la curvatura local del perfil de concentración.

Cuando la distribución presenta un máximo local, su segunda derivada es negativa, siempre que el perfil sea suficientemente regular. La ecuación predice entonces que la concentración en ese punto tiende a disminuir.

Por el contrario, en un mínimo local, donde la segunda derivada es positiva, la concentración tiende a aumentar.

DisminuyeAumentaMáximo localMínimo localccxx
Figura 3. Interpretación local de la ecuación de difusión. En un máximo de concentración con curvatura negativa, la concentración disminuye; en un mínimo con curvatura positiva, aumenta. Los perfiles son esquemáticos y no representan resultados de una simulación.

De esta manera, la ecuación describe la tendencia de la difusión a suavizar las diferencias espaciales de concentración.

4.4. ¿Qué información necesitamos para resolver la ecuación?

La ecuación de difusión establece cómo evoluciona la concentración, pero no determina por sí sola una solución única.

Para predecir el comportamiento de un sistema concreto, también necesitamos especificar su estado inicial y las condiciones que se cumplen en sus fronteras.

La condición inicial establece cómo se encuentra distribuida la sustancia en el instante en que comienza nuestro estudio:

\[c(x,0)=c_0(x),\]

donde \(c_0(x)\) representa el perfil inicial de concentración.

Las condiciones de frontera describen, por su parte, lo que sucede en los extremos del dominio.

Por ejemplo, si el recipiente está cerrado y sus paredes son impermeables, no puede existir flujo de materia a través de ellas. Para un recipiente de longitud \(L\), esto significa que

\[J(0,t)=J(L,t)=0.\]

Si el coeficiente de difusión es positivo, la primera ley de Fick permite expresar estas condiciones como

\[\left.\frac{\partial c}{\partial x}\right|_{x=0}=\left.\frac{\partial c}{\partial x}\right|_{x=L}=0.\]

Estas condiciones tienen una consecuencia importante: la cantidad total de sustancia contenida en el recipiente permanece constante.

En efecto, integrando la ecuación de continuidad a lo largo del dominio y multiplicando por el área transversal, obtenemos

\[\frac{d}{dt}\left(A\int_0^L c(x,t)\,dx\right)=A[J(0,t)-J(L,t)].\]

Como ambos flujos de frontera son nulos,

\[\boxed{\frac{dn_{\mathrm{total}}}{dt}=0}\]

La cantidad total de sustancia se conserva, aunque su distribución espacial cambie continuamente.

Así, hemos llegado a un modelo que combina una ley de transporte, un principio de conservación y las condiciones particulares del sistema que deseamos estudiar.

En las siguientes secciones utilizaremos este modelo para investigar cómo evoluciona una distribución inicial de concentraciones y cómo podemos reproducir esa evolución mediante una simulación numérica.

5. Del movimiento microscópico a la caminata aleatoria

La ecuación de difusión describe la evolución de la concentración sin necesidad de conocer la trayectoria de cada partícula. Esta descripción es particularmente útil cuando estudiamos sistemas formados por un número enorme de partículas, como las moléculas de una sustancia disuelta en agua.

Sin embargo, podemos preguntarnos qué relación existe entre esta descripción macroscópica y el comportamiento de las partículas individuales.

A escala microscópica, las moléculas de un líquido experimentan un movimiento térmico continuo y mantienen interacciones con las partículas que las rodean. Como consecuencia, una partícula disuelta cambia repetidamente la dirección y la magnitud de su velocidad.

En lugar de intentar reconstruir todos los detalles de estas interacciones, construiremos un modelo estadístico simplificado: la caminata aleatoria.

Este modelo no pretende reproducir exactamente la dinámica molecular. Su propósito es representar las propiedades estadísticas fundamentales del movimiento mediante una sucesión de desplazamientos determinados por reglas probabilísticas.

5.1. Una caminata aleatoria en dos dimensiones

Imaginemos dos partículas que se encuentran inicialmente en un mismo punto de un plano. Ambas obedecen las mismas reglas de movimiento, pero sus desplazamientos son independientes.

Durante cada intervalo de tiempo, las partículas pueden desplazarse una distancia fija en una de cuatro direcciones: hacia arriba, hacia abajo, hacia la derecha o hacia la izquierda. Supondremos que las cuatro posibilidades tienen la misma probabilidad.

Cada partícula realiza así una sucesión de desplazamientos cuya dirección se determina aleatoriamente.

Aunque las dos comiencen en el mismo punto y estén sometidas a las mismas reglas estadísticas, generalmente describirán trayectorias diferentes. Es posible que en algunos momentos se aproximen y que en otros se alejen considerablemente.

No existe una dirección privilegiada para ninguna de ellas. Sin embargo, esto no significa que las partículas deban permanecer cerca de su posición inicial. A medida que aumenta el número de pasos, también aumenta, en promedio, el cuadrado de su desplazamiento respecto del punto de partida.

SIMULACIÓN INTERACTIVA Dos caminatas aleatorias en dos dimensiones
Cómo usar esta simulación
  • Pulsa Iniciar para ejecutar la caminata automáticamente.
  • Usa Un paso para observar el proceso lentamente.
  • Pulsa Reiniciar para comenzar nuevamente desde el origen.
  • Observa que ambas partículas obedecen las mismas reglas, pero producen trayectorias distintas.
0 / 200 pasos
● Partícula A● Partícula BOrigen (0, 0)

Las trayectorias son realizaciones independientes del modelo; cada paso tiene una longitud fija y una de cuatro direcciones igualmente probables. La animación se detiene al alcanzar 200 pasos.

Esta sencilla representación nos permite identificar una propiedad importante: las reglas probabilísticas pueden ser las mismas para todas las partículas, aunque sus trayectorias individuales sean diferentes.

Para describir el comportamiento colectivo no necesitamos predecir cada trayectoria. Podemos estudiar magnitudes estadísticas, como la probabilidad de encontrar una partícula en una determinada región o el desplazamiento cuadrático medio de un conjunto de partículas.

5.2. El modelo unidimensional

Para desarrollar las ecuaciones utilizaremos una versión más sencilla de la caminata aleatoria, en la que el movimiento está restringido a una línea recta.

Sea una partícula cuya posición inicial es \(X_0=0\). En cada intervalo de tiempo \(\Delta t\), realiza un desplazamiento de longitud \(\Delta x\) hacia la derecha o hacia la izquierda, con la misma probabilidad:

\[P(\text{derecha})=\frac12,\qquad P(\text{izquierda})=\frac12.\]

Supondremos que los desplazamientos sucesivos son independientes y que las probabilidades no dependen de la posición ni del instante considerado.

Podemos representar el desplazamiento realizado durante el paso \(i\) mediante una variable aleatoria \(\xi_i\), cuyos dos valores posibles son

\[\xi_i=\begin{cases}+\Delta x,&\text{con probabilidad }1/2,\\ -\Delta x,&\text{con probabilidad }1/2.\end{cases}\]

Después de \(N\) pasos, la posición de la partícula será

\[X_N=X_0+\sum_{i=1}^{N}\xi_i.\]

Esta expresión muestra que la posición final depende de la suma de todos los desplazamientos realizados.

Aunque no conocemos de antemano el valor de cada \(\xi_i\), sí conocemos sus propiedades estadísticas. Por ejemplo, su valor esperado es

\[\langle\xi_i\rangle=\frac12\Delta x-\frac12\Delta x=0.\]

Por tanto, el desplazamiento esperado después de \(N\) pasos también es nulo:

\[\langle X_N-X_0\rangle=0.\]

Este resultado expresa que el modelo no favorece ninguna de las dos direcciones.

Sin embargo, un desplazamiento medio nulo no significa que las partículas permanezcan inmóviles ni que regresen necesariamente al punto de partida. Una partícula puede encontrarse muy lejos de su posición inicial, mientras que otra puede haber regresado a ella después del mismo número de pasos.

Para cuantificar esta dispersión necesitaremos estudiar el cuadrado del desplazamiento, como haremos más adelante.

5.3. De las trayectorias individuales a la distribución de probabilidad

Hemos descrito cómo cambia la posición de una partícula, pero nuestro objetivo es construir una descripción estadística que no dependa de una trayectoria particular.

Para ello, imaginemos un conjunto muy grande de partículas independientes que realizan caminatas aleatorias bajo las mismas condiciones.

Aunque inicialmente todas se encuentren en un mismo punto, después de varios pasos estarán distribuidas entre diferentes posiciones.

Denotaremos mediante \(P_j(N)\) la probabilidad de encontrar una partícula en la posición \(x_j=j\Delta x\) después de \(N\) pasos.

Esta distribución contiene información sobre el comportamiento del conjunto. Nos indica qué posiciones son más probables y cómo cambia la dispersión a medida que aumenta el número de pasos.

El siguiente problema consiste en determinar cómo evoluciona matemáticamente esa distribución. Para resolverlo, utilizaremos las reglas probabilísticas que acabamos de establecer.

6. De la caminata aleatoria a la ecuación de difusión

En las secciones anteriores obtuvimos la ecuación de difusión a partir de la primera ley de Fick y del principio de conservación de la materia. Este procedimiento utiliza magnitudes macroscópicas, como la concentración y el flujo, sin necesidad de describir el comportamiento individual de las partículas.

Ahora seguiremos un camino diferente. Partiremos de las reglas probabilísticas de una caminata aleatoria y estudiaremos cómo evoluciona la distribución de posiciones de un conjunto de partículas.

Nuestro objetivo es demostrar que, bajo determinadas condiciones, ambos procedimientos conducen a la misma ecuación matemática.

6.1. La ecuación maestra de la caminata aleatoria

Retomemos el modelo unidimensional de la sección anterior. Una partícula realiza desplazamientos de longitud \(\Delta x\), hacia la derecha o hacia la izquierda, durante intervalos de tiempo \(\Delta t\). Ambos desplazamientos tienen una probabilidad de \(1/2\) y son independientes de los pasos anteriores.

Denotaremos mediante \(P_j(N)\) la probabilidad de encontrar la partícula en la posición

\[x_j=j\Delta x\]

después de \(N\) pasos.

Para que la partícula llegue a la posición \(x_j\) en el paso \(N+1\), solo existen dos posibilidades: que se encontrara en \(x_{j-1}\) y se desplazara hacia la derecha, o que se encontrara en \(x_{j+1}\) y se desplazara hacia la izquierda.

1/21/2xⱼ₋₁xⱼxⱼ₊₁Paso NPaso N+1Paso N
Figura 4. Una partícula puede llegar a \(x_j\) desde cualquiera de sus posiciones vecinas. Las probabilidades corresponden a los desplazamientos condicionales, no a las probabilidades de ocupación.

Como ambas posibilidades son mutuamente excluyentes, la probabilidad de encontrarse en \(x_j\) después de \(N+1\) pasos es

\[\boxed{P_j(N+1)=\frac12P_{j-1}(N)+\frac12P_{j+1}(N)}\]

Esta relación recibe el nombre de ecuación maestra de nuestra caminata aleatoria.

Es importante observar que no nos indica cuál será la trayectoria de una partícula concreta. En cambio, determina cómo evoluciona la distribución de probabilidad a partir de la distribución del paso anterior.

Si conocemos todas las probabilidades \(P_j(N)\), podemos calcular las correspondientes al paso siguiente y repetir el procedimiento tantas veces como sea necesario.

La evolución de la distribución es, por tanto, determinista, aunque las trayectorias individuales sean aleatorias.

6.2. De las diferencias finitas a las derivadas

Nuestro siguiente objetivo consiste en obtener una ecuación que describa la evolución de la distribución mediante variables espaciales y temporales continuas.

Para ello, identificaremos el tiempo transcurrido después de \(N\) pasos como

\[t=N\Delta t.\]

Introduciremos también una función continua \(p(x,t)\), que representa la densidad de probabilidad de encontrar una partícula cerca de la posición \(x\) en el instante \(t\).

En el modelo discreto, la probabilidad de ocupar un sitio de la red es adimensional. En la descripción continua, en cambio, \(p(x,t)\) es una densidad de probabilidad, por lo que tiene unidades de longitud inversa:

\[[p]=\mathrm{m^{-1}}.\]

La probabilidad de encontrar la partícula en un intervalo pequeño de longitud \(\Delta x\) puede aproximarse mediante

\[P_j(N)\simeq p(x_j,t)\Delta x.\]

Sustituyendo esta relación en la ecuación maestra y cancelando el factor común \(\Delta x\), obtenemos

\[p(x,t+\Delta t)=\frac12p(x-\Delta x,t)+\frac12p(x+\Delta x,t).\]

Reordenamos la expresión para separar la variación temporal:

\[p(x,t+\Delta t)-p(x,t)=\frac12\left[p(x-\Delta x,t)-2p(x,t)+p(x+\Delta x,t)\right].\]

Finalmente, dividimos ambos miembros entre \(\Delta t\):

\[\frac{p(x,t+\Delta t)-p(x,t)}{\Delta t}=\frac{(\Delta x)^2}{2\Delta t}\frac{p(x+\Delta x,t)-2p(x,t)+p(x-\Delta x,t)}{(\Delta x)^2}.\]

Esta expresión contiene dos cocientes de diferencias que podemos relacionar con derivadas.

Si la función \(p(x,t)\) es suficientemente regular, cuando \(\Delta t\) es pequeño, el cociente del miembro izquierdo aproxima su derivada temporal:

\[\frac{p(x,t+\Delta t)-p(x,t)}{\Delta t}\longrightarrow\frac{\partial p}{\partial t}.\]

De manera semejante, el cociente de diferencias del miembro derecho aproxima la segunda derivada espacial:

\[\frac{p(x+\Delta x,t)-2p(x,t)+p(x-\Delta x,t)}{(\Delta x)^2}\longrightarrow\frac{\partial^2p}{\partial x^2}.\]

Por tanto, si existe un límite continuo apropiado, podemos obtener una ecuación diferencial que describe la evolución de la densidad de probabilidad.

6.3. El límite difusivo

Hay un detalle importante que debemos examinar antes de tomar el límite: el factor que multiplica a la segunda derivada espacial es

\[\frac{(\Delta x)^2}{2\Delta t}.\]

Si hacemos que \(\Delta x\) y \(\Delta t\) tiendan a cero de manera independiente, este cociente puede tender a cero, crecer indefinidamente o aproximarse a un valor finito.

Para obtener una descripción difusiva no trivial, consideraremos un límite en el que ambas cantidades tienden a cero, pero su relación permanece constante:

\[\boxed{D=\frac{(\Delta x)^2}{2\Delta t}}\]

Esta relación define el coeficiente de difusión del modelo unidimensional.

El límite que estamos construyendo no consiste, por tanto, en reducir arbitrariamente el tamaño de los pasos y los intervalos de tiempo. Debemos hacerlo conservando la escala física que determina la dispersión estadística de las partículas.

Bajo esta condición, y suponiendo que las distribuciones discretas convergen apropiadamente a una densidad suficientemente regular, obtenemos

\[\boxed{\frac{\partial p}{\partial t}=D\frac{\partial^2p}{\partial x^2}}\]

Hemos llegado a una ecuación de difusión para la densidad de probabilidad.

Este resultado establece una conexión entre dos niveles de descripción: las reglas estocásticas que gobiernan los desplazamientos individuales y la ecuación diferencial determinista que describe la evolución de la distribución.

Conviene precisar que el límite continuo no es una descripción exacta de todos los detalles microscópicos. Se trata de un modelo efectivo que reproduce determinadas propiedades estadísticas del movimiento a escalas espaciales y temporales apropiadas.

6.4. De la probabilidad a la concentración

La ecuación que acabamos de obtener describe la densidad de probabilidad de encontrar una partícula en una posición determinada. Sin embargo, en las primeras secciones del laboratorio utilizamos la concentración molar \(c(x,t)\).

¿En qué condiciones podemos relacionar ambas descripciones?

Consideremos un conjunto de partículas independientes e idénticas, sometidas a las mismas reglas de movimiento. Supondremos que sus posiciones están descritas por una misma densidad de probabilidad \(p(x,t)\) y que no existen interacciones que modifiquen significativamente el coeficiente de difusión.

Sea \(n_{\mathrm{total}}\) la cantidad total de sustancia, expresada en moles, contenida en un recipiente de sección transversal constante \(A\).

Como \(p(x,t)\) es una densidad de probabilidad normalizada, la cantidad esperada de sustancia que se encuentra en un pequeño intervalo de longitud \(dx\) será

\[dn=n_{\mathrm{total}}p(x,t)\,dx.\]

Por otra parte, el volumen correspondiente a ese intervalo es \(A\,dx\), por lo que

\[dn=c(x,t)A\,dx.\]

Al comparar ambas expresiones encontramos que

\[c(x,t)=\frac{n_{\mathrm{total}}}{A}p(x,t).\]

Si la cantidad total de sustancia y el área transversal permanecen constantes, el factor de proporcionalidad no depende de la posición ni del tiempo.

En consecuencia,

\[\frac{\partial c}{\partial t}=\frac{n_{\mathrm{total}}}{A}\frac{\partial p}{\partial t},\]

y

\[\frac{\partial^2c}{\partial x^2}=\frac{n_{\mathrm{total}}}{A}\frac{\partial^2p}{\partial x^2}.\]

Sustituyendo la ecuación de evolución de la densidad de probabilidad, obtenemos

\[\boxed{\frac{\partial c}{\partial t}=D\frac{\partial^2c}{\partial x^2}}\]

Esta es precisamente la ecuación de difusión que dedujimos anteriormente mediante la primera ley de Fick y la conservación de la materia.

Hemos recuperado la misma estructura matemática utilizando dos procedimientos diferentes.

El primero parte de magnitudes macroscópicas y de una relación constitutiva entre el flujo y el gradiente de concentración. El segundo parte de un modelo microscópico probabilístico y obtiene una descripción continua de la distribución colectiva.

La coincidencia entre ambos resultados muestra cómo ciertas propiedades macroscópicas pueden comprenderse a partir de modelos estadísticos del movimiento microscópico.

6.5. ¿Qué nos enseña esta conexión?

La deducción anterior ilustra una idea que aparece en numerosos problemas de la física: una descripción macroscópica puede ser determinista incluso cuando el modelo microscópico que la fundamenta es probabilístico.

Esta distinción no significa que las fluctuaciones desaparezcan. En un sistema formado por un número finito de partículas, la distribución observada presenta fluctuaciones alrededor de sus valores esperados.

Sin embargo, cuando el número de partículas es suficientemente grande y sus correlaciones no son demasiado importantes, las fluctuaciones relativas de muchas magnitudes colectivas pueden hacerse pequeñas. En esas condiciones, una descripción continua basada en valores medios puede proporcionar predicciones macroscópicas muy precisas.

También debemos reconocer los límites de nuestro modelo. Hemos supuesto que los pasos son independientes, que no existe una dirección preferente y que sus longitudes e intervalos de tiempo están fijados. Si modificamos estas hipótesis, podemos obtener procesos de transporte con propiedades diferentes.

Por ejemplo, la presencia de fuerzas externas puede introducir una tendencia sistemática al desplazamiento, mientras que las correlaciones temporales entre pasos pueden alterar la manera en que crece la dispersión de las partículas.

Estas posibilidades serán importantes cuando estudiemos el transporte de iones en presencia de campos eléctricos y otros fenómenos relacionados con el funcionamiento de las neuronas.

Por ahora, hemos establecido el resultado fundamental: bajo las hipótesis del modelo y en el límite difusivo, una caminata aleatoria conduce a la misma ecuación de evolución que obtuvimos a partir de la primera ley de Fick.

7. El desplazamiento cuadrático medio y el coeficiente de difusión

En las secciones anteriores hemos construido un modelo de caminata aleatoria y demostrado que, bajo determinadas condiciones, la evolución de su distribución de probabilidad conduce a la ecuación de difusión.

Sin embargo, todavía podemos plantear una pregunta importante: ¿cómo podemos cuantificar la dispersión de las partículas a medida que transcurre el tiempo?

Sabemos que, en una caminata aleatoria simétrica, el desplazamiento medio es nulo. No obstante, las partículas pueden alejarse considerablemente de sus posiciones iniciales.

Para caracterizar este comportamiento necesitamos una magnitud que no dependa del signo del desplazamiento y que permita medir la extensión espacial de la distribución.

Esa magnitud es el desplazamiento cuadrático medio, que desempeña un papel fundamental en el estudio de la difusión y de numerosos procesos estocásticos.

7.1. ¿Por qué el desplazamiento medio no es suficiente?

Retomemos nuestra caminata aleatoria unidimensional. La posición de una partícula después de \(N\) pasos está dada por

\[X_N=X_0+\sum_{i=1}^{N}\xi_i,\]

donde \(\xi_i\) representa el desplazamiento correspondiente al paso \(i\).

Recordemos que cada desplazamiento puede tomar los valores \(+\Delta x\) o \(-\Delta x\), ambos con probabilidad \(1/2\).

Por tanto, su valor esperado es

\[\langle\xi_i\rangle=0.\]

El desplazamiento total respecto de la posición inicial es

\[X_N-X_0=\sum_{i=1}^{N}\xi_i.\]

Al calcular su valor esperado y utilizar la linealidad de la esperanza matemática, obtenemos

\[\langle X_N-X_0\rangle=\sum_{i=1}^{N}\langle\xi_i\rangle=0.\]

Este resultado expresa que, al considerar un conjunto suficientemente grande de caminatas independientes, los desplazamientos hacia la derecha y hacia la izquierda se compensan estadísticamente.

Sin embargo, no proporciona información sobre la distancia que las partículas se han alejado del punto de partida.

Imaginemos, por ejemplo, dos partículas que parten del origen. Después de cierto tiempo, una se encuentra en \(x=10\,\mu\mathrm{m}\) y la otra en \(x=-10\,\mu\mathrm{m}\).

El desplazamiento medio de ambas es cero, aunque ninguna permanece en su posición inicial.

Necesitamos, por tanto, una medida que caracterice la dispersión espacial sin que los desplazamientos de signos opuestos se cancelen.

7.2. El desplazamiento cuadrático medio

Una manera de resolver este problema consiste en elevar al cuadrado el desplazamiento de cada partícula y calcular después su valor esperado.

Definimos así el desplazamiento cuadrático medio, que abreviaremos como DCM:

\[\boxed{\operatorname{DCM}(N)=\left\langle(X_N-X_0)^2\right\rangle}\]

En la literatura científica en inglés, esta magnitud se conoce como mean squared displacement (MSD).

Su unidad en el Sistema Internacional es el metro cuadrado:

\[[\operatorname{DCM}]=\mathrm{m^2}.\]

El desplazamiento cuadrático medio no es una distancia, sino una medida cuadrática de la dispersión respecto de la posición inicial.

Su raíz cuadrada, en cambio, tiene unidades de longitud y proporciona una escala característica del desplazamiento:

\[\ell_{\mathrm{rms}}(N)=\sqrt{\left\langle(X_N-X_0)^2\right\rangle}.\]

En ocasiones, esta cantidad recibe el nombre de desplazamiento cuadrático medio en su forma de raíz, aunque conviene distinguirla del propio DCM para evitar ambigüedades.

7.3. Deducción a partir de la caminata aleatoria

Para calcular el desplazamiento cuadrático medio, partiremos de la expresión que obtuvimos para el desplazamiento total:

\[X_N-X_0=\sum_{i=1}^{N}\xi_i.\]

Elevando al cuadrado ambos miembros,

\[(X_N-X_0)^2=\left(\sum_{i=1}^{N}\xi_i\right)^2.\]

Al desarrollar el cuadrado de la suma, aparecen dos tipos de términos:

\[(X_N-X_0)^2=\sum_{i=1}^{N}\xi_i^2+2\sum_{i<j}\xi_i\xi_j.\]

El primer término contiene los cuadrados de los desplazamientos individuales. El segundo contiene los productos entre desplazamientos correspondientes a pasos diferentes.

Calculamos ahora el valor esperado:

\[\left\langle(X_N-X_0)^2\right\rangle=\sum_{i=1}^{N}\langle\xi_i^2\rangle+2\sum_{i<j}\langle\xi_i\xi_j\rangle.\]

Examinemos por separado ambos términos.

Como cada paso tiene una longitud fija \(\Delta x\), independientemente de su dirección,

\[\xi_i^2=(\Delta x)^2.\]

En consecuencia,

\[\langle\xi_i^2\rangle=(\Delta x)^2.\]

Por otra parte, hemos supuesto que los desplazamientos sucesivos son independientes. Para dos pasos diferentes, esta propiedad nos permite escribir

\[\langle\xi_i\xi_j\rangle=\langle\xi_i\rangle\langle\xi_j\rangle,\qquad i\ne j.\]

Como el valor esperado de cada desplazamiento es cero, todos los términos cruzados desaparecen:

\[\langle\xi_i\xi_j\rangle=0,\qquad i\ne j.\]

De esta manera, únicamente permanecen los \(N\) términos correspondientes a los cuadrados de los desplazamientos individuales:

\[\boxed{\left\langle(X_N-X_0)^2\right\rangle=N(\Delta x)^2}\]

Hemos obtenido un resultado importante: el desplazamiento cuadrático medio crece linealmente con el número de pasos.

Este comportamiento no depende de una trayectoria particular. Es una propiedad estadística del conjunto de caminatas aleatorias que obedecen las hipótesis de nuestro modelo.

7.4. El tiempo y el coeficiente de difusión

En nuestro modelo, cada paso requiere un intervalo de tiempo \(\Delta t\). Por tanto, después de \(N\) pasos, el tiempo transcurrido es

\[t=N\Delta t.\]

Podemos despejar el número de pasos:

\[N=\frac{t}{\Delta t}.\]

Sustituyendo esta expresión en el resultado anterior,

\[\left\langle(X_N-X_0)^2\right\rangle=\frac{(\Delta x)^2}{\Delta t}\,t.\]

Recordemos que, al estudiar el límite continuo de la caminata aleatoria, identificamos el coeficiente de difusión mediante la relación

\[D=\frac{(\Delta x)^2}{2\Delta t}.\]

En consecuencia,

\[\frac{(\Delta x)^2}{\Delta t}=2D.\]

Obtenemos finalmente

\[\boxed{\left\langle[X(t)-X(0)]^2\right\rangle=2Dt}\]

Esta relación caracteriza la difusión normal en una dimensión, dentro de las condiciones establecidas por nuestro modelo.

El resultado proporciona una interpretación física del coeficiente de difusión: \(D\) determina la rapidez con la que crece el desplazamiento cuadrático medio.

Si dos sustancias tienen coeficientes de difusión diferentes, sus distribuciones espaciales se ensancharán a ritmos diferentes, siempre que las demás condiciones sean comparables.

Es importante observar que el desplazamiento característico no crece linealmente con el tiempo. Al tomar la raíz cuadrada,

\[\boxed{\ell_{\mathrm{rms}}(t)=\sqrt{2Dt}}\]

encontramos que la escala espacial característica del desplazamiento crece proporcionalmente a la raíz cuadrada del tiempo.

Por ejemplo, para duplicar esa escala espacial característica, el tiempo necesario debe multiplicarse por cuatro.

Esta dependencia tiene consecuencias importantes para el transporte de sustancias en los sistemas biológicos, donde la difusión puede ser eficaz a distancias microscópicas, pero resultar mucho más lenta al considerar distancias considerablemente mayores.

7.5. Generalización a dos y tres dimensiones

Hasta ahora hemos trabajado con una caminata aleatoria unidimensional. Sin embargo, las partículas de un fluido real pueden desplazarse en las tres dimensiones del espacio.

Podemos generalizar nuestro resultado suponiendo que el proceso es isótropo, es decir, que sus propiedades estadísticas son las mismas en todas las direcciones.

En dos dimensiones, el vector de desplazamiento es

\[\Delta\mathbf r=(\Delta X,\Delta Y).\]

El cuadrado de su magnitud es

\[|\Delta\mathbf r|^2=(\Delta X)^2+(\Delta Y)^2.\]

Por linealidad de la esperanza matemática,

\[\left\langle|\Delta\mathbf r|^2\right\rangle=\langle(\Delta X)^2\rangle+\langle(\Delta Y)^2\rangle.\]

Para una difusión isótropa con el mismo coeficiente \(D\) en ambas direcciones, cada coordenada satisface

\[\langle(\Delta X)^2\rangle=\langle(\Delta Y)^2\rangle=2Dt.\]

Por tanto,

\[\boxed{\left\langle|\Delta\mathbf r|^2\right\rangle=4Dt}\]

En tres dimensiones, el mismo razonamiento conduce a

\[\boxed{\left\langle|\Delta\mathbf r|^2\right\rangle=6Dt}\]

De manera más general, para una difusión normal e isótropa en un espacio de \(d\) dimensiones,

\[\boxed{\left\langle|\Delta\mathbf r|^2\right\rangle=2dDt}\]

Esta relación permite interpretar los resultados de nuestra simulación bidimensional de la sección 5. Aunque las trayectorias de las dos partículas sean diferentes, la predicción estadística para un conjunto numeroso de caminatas independientes es que el desplazamiento cuadrático medio crecerá proporcionalmente al tiempo, con una pendiente igual a \(4D\).

0750150022503000Tiempo (s)DCM (µm²)0123451D: 2Dt2D: 4Dt3D: 6Dt
Figura 5. Crecimiento teórico del desplazamiento cuadrático medio para difusión normal e isótropa en una, dos y tres dimensiones, con \(D=100\,\mu\mathrm{m}^2/\mathrm{s}\). Las rectas no son datos experimentales ni resultados de simulación.

7.6. Una predicción que podemos comprobar

La relación entre el desplazamiento cuadrático medio y el tiempo nos proporciona una predicción cuantitativa que podemos contrastar mediante simulaciones numéricas.

Si generamos un conjunto suficientemente numeroso de caminatas aleatorias independientes y registramos las posiciones de las partículas en distintos instantes, podremos calcular el desplazamiento cuadrático medio del conjunto:

\[\operatorname{DCM}_M(t)=\frac1M\sum_{k=1}^{M}\left|\mathbf r_k(t)-\mathbf r_k(0)\right|^2,\]

donde \(M\) representa el número de partículas y \(\mathbf r_k(t)\) es la posición de la partícula \(k\) en el instante \(t\).

Esta expresión es un estimador del desplazamiento cuadrático medio teórico obtenido a partir de la distribución de probabilidad.

Para un conjunto finito de partículas, los valores calculados presentarán fluctuaciones estadísticas. Por ello, no debemos esperar que los resultados coincidan exactamente con la predicción teórica en todos los instantes.

Sin embargo, al aumentar el número de caminatas independientes, el estimador tiende a aproximarse al valor esperado.

En una simulación bidimensional, esperamos encontrar una relación aproximadamente lineal:

\[\operatorname{DCM}_M(t)\simeq 4Dt.\]

Podemos utilizar esta relación para estimar el coeficiente de difusión a partir de la pendiente de una gráfica del desplazamiento cuadrático medio frente al tiempo.

Así, una magnitud que introdujimos para caracterizar la dispersión estadística se convierte también en una herramienta para determinar un parámetro físico del modelo.

Esta será precisamente una de las preguntas que investigaremos mediante nuestra siguiente simulación numérica.

8. Soluciones analíticas de la ecuación de difusión

Hasta ahora hemos obtenido la ecuación de difusión mediante dos procedimientos diferentes. El primero se basó en la conservación de la materia y en la primera ley de Fick; el segundo partió de un modelo de caminata aleatoria y de su límite continuo. Ambos procedimientos condujeron a la misma ecuación:

\[\boxed{\frac{\partial c}{\partial t}=D\frac{\partial^2c}{\partial x^2}}\]

Esta ecuación establece cómo se relaciona la variación temporal de la concentración con la curvatura de su distribución espacial. Sin embargo, por sí sola no nos indica cuál es la concentración en cada posición e instante. Para obtener esa información debemos resolver la ecuación, incorporando las condiciones iniciales y de frontera que caracterizan el sistema físico.

En esta sección estudiaremos dos soluciones analíticas especialmente importantes: la evolución de una cantidad de sustancia inicialmente localizada y la de una distribución que posee una anchura inicial finita. Ambas nos permitirán comprender cómo se dispersa la materia, cómo se conserva su cantidad total y qué relación existe entre la difusión y el desplazamiento cuadrático medio. Finalmente, examinaremos cómo se modifica el problema cuando el transporte ocurre en un recipiente de longitud finita con paredes impermeables.

8.1. El problema de una cantidad de sustancia inicialmente localizada

Consideremos un medio homogéneo, de sección transversal constante \(A\), que se extiende idealmente desde \(x=-\infty\) hasta \(x=+\infty\). Supondremos que el coeficiente de difusión \(D\) es constante y que no existen corrientes macroscópicas, reacciones químicas ni otros mecanismos de transporte.

En el instante inicial, toda la cantidad de sustancia, \(n_0\), se encuentra concentrada en una región que idealizaremos como un punto situado en el origen. Matemáticamente, esta condición inicial se representa mediante

\[c(x,0)=\frac{n_0}{A}\delta(x),\]

donde \(\delta(x)\) es la delta de Dirac. La delta de Dirac no es una función ordinaria. Es una distribución matemática que permite representar una cantidad finita concentrada en un punto. Su propiedad fundamental, en este contexto, es

\[\int_{-\infty}^{+\infty}\delta(x)\,dx=1.\]

Por tanto, la cantidad inicial de sustancia es

\[A\int_{-\infty}^{+\infty}c(x,0)\,dx=n_0.\]

La condición inicial es una idealización: en un experimento real, la sustancia siempre ocupa una región de volumen finito. Sin embargo, esta representación resulta especialmente útil cuando queremos estudiar la respuesta del sistema a una distribución inicialmente muy localizada.

Nuestro problema consiste en encontrar una función \(c(x,t)\) que satisfaga simultáneamente la ecuación de difusión y esa condición inicial, suponiendo además que la concentración y su gradiente desaparecen suficientemente rápido al alejarnos hacia el infinito.

8.2. Resolución mediante la transformada de Fourier

La dificultad de nuestra ecuación es que contiene derivadas respecto de dos variables independientes: la posición y el tiempo. Podemos simplificar el problema mediante una herramienta matemática que permite representar una distribución espacial como una superposición de componentes oscilatorias: la transformada de Fourier.

Para una función \(c(x,t)\), definiremos su transformada espacial mediante

\[\widehat c(k,t)=\int_{-\infty}^{+\infty}c(x,t)e^{-ikx}\,dx,\]

donde \(k\) representa el número de onda, cuyas unidades son \(\mathrm{m^{-1}}\). La transformada inversa permite reconstruir la distribución original:

\[c(x,t)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}\widehat c(k,t)e^{ikx}\,dk.\]

Estas dos expresiones establecen una correspondencia entre la distribución espacial de la concentración y su representación en términos de números de onda. La ventaja fundamental es que, bajo las condiciones de regularidad y decaimiento apropiadas, la transformada convierte la segunda derivada espacial en una multiplicación por \(-k^2\):

\[\mathcal F\!\left[\frac{\partial^2c}{\partial x^2}\right]=-k^2\widehat c(k,t).\]

Esta propiedad puede obtenerse integrando por partes dos veces, utilizando que los términos de frontera desaparecen. Al aplicar la transformada de Fourier a la ecuación de difusión, obtenemos

\[\frac{\partial\widehat c}{\partial t}=-Dk^2\widehat c.\]

Hemos transformado una ecuación diferencial parcial en una ecuación diferencial ordinaria respecto del tiempo, para cada número de onda \(k\). Su solución es

\[\widehat c(k,t)=\widehat c(k,0)e^{-Dk^2t}.\]

Podemos interpretar este resultado antes de continuar con la deducción. Cada componente espacial evoluciona de manera independiente y su amplitud disminuye exponencialmente con el tiempo. Además, las componentes con números de onda de mayor magnitud decaen más rápidamente que las correspondientes a números de onda pequeños.

Esto significa que las variaciones espaciales que ocurren a escalas pequeñas desaparecen más rápidamente que las variaciones de gran escala. Es una expresión matemática de la tendencia de la difusión a suavizar los perfiles de concentración.

Incorporación de la condición inicial

Para nuestra distribución inicialmente localizada, \(c(x,0)=\frac{n_0}{A}\delta(x)\), la transformada de Fourier es

\[\widehat c(k,0)=\frac{n_0}{A}\int_{-\infty}^{+\infty}\delta(x)e^{-ikx}\,dx.\]

Utilizando la propiedad de la delta de Dirac, \(\widehat c(k,0)=\frac{n_0}{A}\). Por tanto,

\[\widehat c(k,t)=\frac{n_0}{A}e^{-Dk^2t}.\]

Para regresar al espacio de posiciones, aplicamos la transformada inversa:

\[c(x,t)=\frac{n_0}{2\pi A}\int_{-\infty}^{+\infty}e^{-Dk^2t}e^{ikx}\,dk.\]

La integral que aparece en esta expresión es una integral gaussiana. Para \(D>0\) y \(t>0\), su resultado es

\[\int_{-\infty}^{+\infty}e^{-Dtk^2+ikx}\,dk=\sqrt{\frac{\pi}{Dt}}\exp\!\left(-\frac{x^2}{4Dt}\right).\]

Sustituyendo esta identidad y simplificando, obtenemos finalmente

\[\boxed{c(x,t)=\frac{n_0}{A\sqrt{4\pi Dt}}\exp\!\left(-\frac{x^2}{4Dt}\right)}\]

Esta expresión es la solución fundamental de la ecuación de difusión unidimensional, multiplicada por la cantidad inicial de sustancia por unidad de área. Es válida para tiempos positivos y satisface la condición inicial localizada en el sentido de las distribuciones: conforme \(t\) tiende a cero por valores positivos, la solución converge hacia la delta de Dirac con el factor correspondiente.

8.3. Interpretación física de la solución gaussiana

La solución que acabamos de obtener contiene información sobre diferentes aspectos del proceso de difusión. Podemos escribirla como

\[c(x,t)=\frac{n_0}{A\sqrt{4\pi Dt}}e^{-x^2/(4Dt)}.\]

Su forma espacial es gaussiana y presenta un máximo en el origen. A medida que transcurre el tiempo, la distribución se ensancha y su altura disminuye.

Posición xConcentración ct₁t₂t₃t₁ < t₂ < t₃
Figura 6. Evolución esquemática de una distribución gaussiana por difusión. Conforme aumenta el tiempo, disminuye la concentración máxima y aumenta la anchura del perfil. El área bajo las tres curvas permanece constante.

Visualización de la evolución temporal

La figura anterior resume la evolución de una distribución gaussiana en tres instantes. La siguiente animación, generada mediante Python a partir de la solución analítica, permite observar de forma continua cómo el perfil se ensancha y se aplana. No representa una simulación de partículas individuales, sino la evolución de la concentración predicha por la ecuación de difusión. La normalización y el dominio espacial mostrado se indican en el programa descargable.

Animación: ensanchamiento y aplanamiento temporal de una distribución gaussiana
Figura 7. Animación de la solución gaussiana para \(D=100\,\mu\mathrm{m}^2/\mathrm{s}\). El perfil se ensancha mientras disminuye su máximo; su integral en el dominio infinito permanece constante.

Descargar la animación GIF ·

Conservación de la materia

La primera propiedad que debemos comprobar es que la solución respeta el principio de conservación. La cantidad total de sustancia contenida en el sistema es

\[n(t)=A\int_{-\infty}^{+\infty}c(x,t)\,dx.\]

Sustituyendo la solución gaussiana,

\[n(t)=\frac{n_0}{\sqrt{4\pi Dt}}\int_{-\infty}^{+\infty}e^{-x^2/(4Dt)}\,dx.\]

La integral gaussiana tiene el valor

\[\int_{-\infty}^{+\infty}e^{-x^2/(4Dt)}\,dx=\sqrt{4\pi Dt}.\]

Los factores se cancelan y obtenemos \(\boxed{n(t)=n_0}\). La cantidad total de sustancia permanece constante, aunque su concentración disminuya en unas regiones y aumente en otras. Esta propiedad tiene una interpretación geométrica: el área bajo la curva de concentración, multiplicada por la sección transversal del medio, es independiente del tiempo.

Evolución de la concentración máxima

La concentración máxima se encuentra en el origen. Al evaluar la solución en \(x=0\), obtenemos

\[c(0,t)=\frac{n_0}{A\sqrt{4\pi Dt}}.\]

Por tanto, \(\boxed{c_{\mathrm{max}}(t)\propto t^{-1/2}}\). La concentración máxima disminuye a medida que la sustancia se dispersa. Esta disminución no se debe a la desaparición de partículas, sino a su distribución en una región espacial progresivamente más extensa.

Anchura de la distribución y desplazamiento cuadrático medio

Podemos establecer una conexión directa con el resultado que obtuvimos en la sección 7. Si normalizamos la concentración respecto de la cantidad total de sustancia, obtenemos la densidad de probabilidad

\[p(x,t)=\frac{A}{n_0}c(x,t)=\frac{1}{\sqrt{4\pi Dt}}\exp\!\left(-\frac{x^2}{4Dt}\right).\]

Esta es una distribución gaussiana de media cero y varianza \(2Dt\). Su segundo momento es

\[\langle X^2(t)\rangle=\int_{-\infty}^{+\infty}x^2p(x,t)\,dx.\]

Al evaluar la integral encontramos \(\boxed{\langle X^2(t)\rangle=2Dt}\). Hemos recuperado exactamente la relación que dedujimos a partir de la caminata aleatoria. Este resultado establece una correspondencia entre la descripción estadística microscópica y la solución analítica de la ecuación de difusión. Ambos procedimientos predicen que la dispersión cuadrática media aumenta linealmente con el tiempo.

La desviación estándar de la distribución es \(\boxed{\sigma(t)=\sqrt{2Dt}}\). Así, el ensanchamiento de la distribución gaussiana tiene la misma dependencia temporal que la escala característica del desplazamiento de las partículas.

8.4. Evolución de una distribución gaussiana de anchura inicial finita

La solución fundamental describe una cantidad de sustancia inicialmente concentrada en un punto. Sin embargo, en un sistema físico real es más habitual que la distribución inicial ocupe una región de anchura finita.

Consideremos una distribución gaussiana centrada en \(x_0\), cuya anchura inicial está caracterizada por el parámetro \(\sigma_0\):

\[c(x,0)=\frac{n_0}{A\sqrt{2\pi\sigma_0^2}}\exp\!\left[-\frac{(x-x_0)^2}{2\sigma_0^2}\right].\]

El parámetro \(\sigma_0\), expresado en unidades de longitud, es la desviación estándar inicial de la distribución. La cantidad total de sustancia es \(n_0\), independientemente de la anchura inicial.

La solución mediante transformadas de Fourier

Podemos utilizar el procedimiento desarrollado anteriormente para resolver este nuevo problema. La transformada de Fourier de la condición inicial es

\[\widehat c(k,0)=\frac{n_0}{A}\exp\!\left(-ikx_0-\frac{\sigma_0^2k^2}{2}\right).\]

Sabemos que cada componente espacial evoluciona según el factor \(e^{-Dk^2t}\). Por tanto,

\[\widehat c(k,t)=\frac{n_0}{A}\exp\!\left[-ikx_0-\left(\frac{\sigma_0^2}{2}+Dt\right)k^2\right].\]

Observemos que podemos escribir el coeficiente cuadrático del exponente como \(\frac{\sigma_0^2}{2}+Dt=\frac{\sigma_0^2+2Dt}{2}\). Al aplicar la transformada inversa obtenemos

\[\boxed{c(x,t)=\frac{n_0}{A\sqrt{2\pi(\sigma_0^2+2Dt)}}\exp\!\left[-\frac{(x-x_0)^2}{2(\sigma_0^2+2Dt)}\right]}\]

El resultado tiene una propiedad notable: la distribución conserva su forma gaussiana durante toda la evolución, pero su anchura aumenta con el tiempo. La varianza está dada por \(\boxed{\sigma^2(t)=\sigma_0^2+2Dt}\). El centro de la distribución permanece en \(x_0\), mientras que su varianza aumenta linealmente.

Si hacemos que la anchura inicial tienda a cero, recuperamos la solución fundamental para una cantidad de sustancia inicialmente localizada en \(x_0\).

Una consecuencia importante

Podemos despejar el tiempo necesario para que la distribución alcance una anchura determinada:

\[t=\frac{\sigma^2(t)-\sigma_0^2}{2D}.\]

Esta expresión permite estimar la escala temporal de la difusión cuando conocemos el coeficiente \(D\) y las dimensiones espaciales del problema. La dependencia cuadrática respecto de la longitud explica por qué el transporte difusivo resulta particularmente eficaz a escalas microscópicas y por qué los tiempos característicos aumentan considerablemente cuando crecen las distancias.

8.5. Difusión en un recipiente con paredes impermeables

Hasta ahora hemos considerado un medio de extensión infinita. Esta aproximación es útil para obtener soluciones analíticas sencillas, pero no siempre representa las condiciones de un sistema físico real.

Consideremos ahora un recipiente de longitud \(L\), cuyos extremos se encuentran en \(x=0\) y \(x=L\). Supondremos que las paredes son impermeables, por lo que no puede existir flujo de materia a través de ellas. Nuestro problema está definido por la ecuación de difusión,

\[\frac{\partial c}{\partial t}=D\frac{\partial^2c}{\partial x^2},\]

y por las condiciones de frontera

\[\left.\frac{\partial c}{\partial x}\right|_{x=0}=\left.\frac{\partial c}{\partial x}\right|_{x=L}=0.\]

Estas condiciones, conocidas como condiciones de Neumann homogéneas, expresan que el flujo difusivo es nulo en ambos extremos. A diferencia de lo que ocurre en el medio infinito, la sustancia no puede dispersarse indefinidamente. Conforme transcurre el tiempo, su distribución tiende a hacerse uniforme dentro del recipiente.

Solución mediante separación de variables

Podemos buscar soluciones de la forma \(c(x,t)=X(x)T(t)\), donde \(X(x)\) describe la dependencia espacial y \(T(t)\) la dependencia temporal. Al sustituir esta expresión en la ecuación de difusión, obtenemos

\[X(x)\frac{dT}{dt}=DT(t)\frac{d^2X}{dx^2}.\]

Dividiendo entre \(DX(x)T(t)\), encontramos

\[\frac{1}{DT}\frac{dT}{dt}=\frac{1}{X}\frac{d^2X}{dx^2}.\]

El miembro izquierdo depende únicamente del tiempo, mientras que el derecho depende únicamente de la posición. Para que la igualdad se cumpla en todos los puntos e instantes, ambos miembros deben ser iguales a una constante. Escribiremos esa constante como \(-\lambda\), con lo que obtenemos dos ecuaciones diferenciales ordinarias:

\[\frac{d^2X}{dx^2}+\lambda X=0,\qquad\frac{dT}{dt}+D\lambda T=0.\]

Las condiciones de frontera exigen que \(X'(0)=X'(L)=0\). Estas condiciones seleccionan las soluciones espaciales

\[X_m(x)=\cos\!\left(\frac{m\pi x}{L}\right),\qquad \lambda_m=\left(\frac{m\pi}{L}\right)^2,\quad m=0,1,2,\ldots\]

Para cada uno de estos valores, la ecuación temporal tiene como solución

\[T_m(t)=\exp\!\left[-D\left(\frac{m\pi}{L}\right)^2t\right].\]

Como la ecuación de difusión es lineal, podemos construir una solución general mediante la superposición de estos modos:

\[\boxed{c(x,t)=\overline c+\sum_{m=1}^{\infty}a_m\cos\!\left(\frac{m\pi x}{L}\right)\exp\!\left[-D\left(\frac{m\pi}{L}\right)^2t\right]}\]

Aquí, \(\overline c\) es la concentración media inicial,

\[\overline c=\frac{1}{L}\int_0^Lc(x,0)\,dx,\]

y los coeficientes \(a_m\) dependen de la distribución inicial:

\[a_m=\frac{2}{L}\int_0^Lc(x,0)\cos\!\left(\frac{m\pi x}{L}\right)\,dx.\]

Cada término de la serie representa un modo espacial cuya amplitud disminuye exponencialmente con el tiempo. El modo \(m=0\), que corresponde a la concentración media, permanece constante. Los modos con \(m\geq1\), responsables de las variaciones espaciales de la concentración, decaen progresivamente.

En consecuencia, \(\boxed{\lim_{t\to\infty}c(x,t)=\overline c}\). El sistema se aproxima a una distribución uniforme, determinada por la cantidad total de sustancia y el volumen del recipiente.

La escala temporal de relajación

La rapidez con la que desaparece cada modo está caracterizada por un tiempo

\[\tau_m=\frac{L^2}{Dm^2\pi^2},\qquad m\geq1.\]

Los modos de mayor número de onda decaen más rápidamente. A tiempos suficientemente largos, el comportamiento de una distribución inicial no uniforme está dominado, por lo general, por el modo de menor número de onda cuyo coeficiente sea distinto de cero. Cuando el coeficiente \(a_1\) no es nulo, el tiempo característico más largo es \(\boxed{\tau_1=\frac{L^2}{\pi^2D}}\).

Esta expresión recupera una característica fundamental del transporte difusivo: el tiempo necesario para suavizar una distribución aumenta con el cuadrado de la longitud característica del sistema. También permite distinguir entre dos regímenes. Cuando los tiempos son suficientemente cortos para que la sustancia no haya alcanzado las paredes, la solución de un medio infinito puede proporcionar una buena aproximación local. A tiempos más largos, las condiciones de frontera se vuelven esenciales para describir correctamente la evolución.

8.6. Conexión con nuestras simulaciones numéricas

Las soluciones que hemos obtenido proporcionan referencias analíticas para interpretar y comprobar nuestros experimentos computacionales. La solución fundamental gaussiana describe la evolución continua de una distribución inicialmente localizada. Será especialmente útil cuando comparemos la ecuación de difusión con el comportamiento colectivo de las caminatas aleatorias.

La solución con anchura inicial finita permite estudiar cómo evoluciona una distribución gaussiana y ofrece una predicción sencilla para su varianza: \(\sigma^2(t)=\sigma_0^2+2Dt\).

Finalmente, la solución mediante una serie de modos espaciales describe la difusión en un recipiente con extremos impermeables. Esta es la geometría que utilizamos en nuestro programa de volúmenes finitos.

El programa numérico emplea una distribución inicial gaussiana superpuesta a una concentración de fondo. Dicha distribución no coincide, en general, con un único modo espacial de la solución analítica. Sin embargo, podemos representarla mediante una serie de cosenos, calcular su evolución y comparar los resultados con la solución numérica.

Esta comparación permitirá examinar si el programa reproduce correctamente la dinámica del modelo, respeta las condiciones de frontera y conserva la cantidad total de sustancia.

Debemos distinguir, no obstante, entre una coincidencia con la solución matemática y una validación experimental del modelo físico. Comparar dos procedimientos matemáticos para resolver el mismo problema permite verificar la implementación numérica, pero no demuestra por sí solo que el modelo represente adecuadamente un sistema biológico real.

Las soluciones analíticas también nos permiten comprender los límites de nuestras aproximaciones. La geometría del dominio, la distribución inicial y las propiedades del medio determinan qué solución debemos utilizar y hasta qué punto sus predicciones son aplicables.

Hemos completado así el recorrido que va desde una ecuación diferencial hasta sus soluciones analíticas y su interpretación física. Disponemos ahora de las herramientas necesarias para estudiar la difusión mediante experimentos computacionales y comparar sus resultados con predicciones cuantitativas.

9. Simulaciones numéricas: de las caminatas aleatorias a la difusión

En las secciones anteriores hemos estudiado la difusión desde diferentes perspectivas. Partimos del movimiento microscópico de las partículas, construimos un modelo de caminata aleatoria, dedujimos la ecuación de difusión y obtuvimos algunas de sus soluciones analíticas.

Ahora queremos examinar hasta qué punto un experimento numérico reproduce las predicciones de esos modelos. Para ello desarrollaremos dos simulaciones complementarias. En la primera representaremos el movimiento de un conjunto de partículas mediante caminatas aleatorias independientes; en la segunda resolveremos directamente la ecuación de difusión mediante un método numérico basado en la conservación de la materia.

La comparación nos permitirá comprender qué información proporciona cada descripción, identificar las fluctuaciones estadísticas y examinar las limitaciones de nuestras aproximaciones.

9.1. Un experimento con mil partículas

Consideremos un conjunto de \(M\) partículas que inicialmente se encuentran en el origen de un eje horizontal. Cada partícula realiza una caminata aleatoria unidimensional y, durante cada intervalo \(\Delta t\), se desplaza una distancia \(\Delta x\) hacia la derecha o hacia la izquierda, con la misma probabilidad.

Las caminatas son independientes y obedecen las mismas reglas estadísticas. Utilizaremos los siguientes parámetros:

ParámetroSímboloValor
Número de partículas\(M\)1 000
Número de pasos\(N\)200
Longitud de cada paso\(\Delta x\)1 µm
Duración de cada paso\(\Delta t\)0.005 s
Posición inicial\(X_0\)0 µm

Estos valores son ilustrativos y no corresponden a mediciones de una sustancia particular. El coeficiente de difusión del modelo es \(D=(\Delta x)^2/(2\Delta t)=100\,\mu\mathrm{m}^2/\mathrm{s}\), y la duración total del experimento es \(T=N\Delta t=1\,\mathrm{s}\).

Construcción del programa

Utilizaremos NumPy para generar los desplazamientos de todas las partículas y calcular sus posiciones. El programa empleará un generador de números aleatorios con semilla fija para que el experimento sea reproducible.

import numpy as np
import matplotlib.pyplot as plt

# Parámetros físicos y numéricos
M = 1000
N = 200
dx = 1.0       # µm
dt = 0.005     # s
D = dx**2 / (2 * dt)

# Generador reproducible
rng = np.random.default_rng(2026)

# Desplazamientos de las partículas
pasos = rng.choice([-dx, dx], size=(M, N))

# Posiciones acumuladas, con la posición inicial
posiciones = np.cumsum(pasos, axis=1)
posiciones = np.column_stack((np.zeros(M), posiciones))
tiempos = np.arange(N + 1) * dt

La matriz pasos contiene una fila por partícula y una columna por desplazamiento. La función np.cumsum calcula las sumas acumuladas y permite obtener la trayectoria completa de cada partícula. La matriz posiciones contiene también la posición inicial, por lo que tiene \(M\) filas y \(N+1\) columnas.

Una semilla fija permite reproducir una realización concreta del experimento. Si modificamos la semilla, obtendremos otras trayectorias, aunque las propiedades estadísticas del modelo permanecerán inalteradas.

fig, ax = plt.subplots(figsize=(9, 5))
for k in range(10):
    ax.plot(tiempos, posiciones[k], linewidth=1)
ax.set_xlabel("Tiempo (s)")
ax.set_ylabel("Posición (µm)")
ax.set_title("Diez caminatas aleatorias")
ax.grid(alpha=0.25)
plt.tight_layout()
plt.show()

Cada curva corresponde a una realización individual del proceso. Algunas partículas se alejarán rápidamente del origen, otras cambiarán frecuentemente de dirección y algunas regresarán a posiciones que ya habían ocupado. El comportamiento de una trayectoria particular no permite determinar por sí solo las propiedades estadísticas del conjunto.

9.2. La distribución de posiciones y la solución gaussiana

En lugar de estudiar las trayectorias individuales, examinemos ahora las posiciones de las mil partículas después de completar los 200 pasos. Construiremos un histograma de sus posiciones finales. En nuestra caminata, después de un número par de pasos, las posiciones accesibles son múltiplos pares de \(\Delta x\), por lo que utilizaremos intervalos de anchura \(2\Delta x\).

finales = posiciones[:, -1]
bordes = np.arange(-N - 1, N + 2, 2) * dx
fig, ax = plt.subplots(figsize=(9, 5))
ax.hist(finales, bins=bordes, density=True,
        alpha=0.6, label="Caminatas aleatorias")
ax.set_xlabel("Posición (µm)")
ax.set_ylabel("Densidad (1/µm)")
ax.set_title("Distribución de posiciones finales")
ax.legend()
plt.tight_layout()
plt.show()

El histograma proporciona una estimación estadística de la distribución de posiciones. Su forma depende de la realización numérica y presenta fluctuaciones porque estamos utilizando un conjunto finito de partículas.

Podemos compararlo con la solución fundamental de la ecuación de difusión obtenida en la sección 8:

x = np.linspace(-60, 60, 600)
p_teorica = (np.exp(-x**2 / (4 * D * T))
             / np.sqrt(4 * np.pi * D * T))
fig, ax = plt.subplots(figsize=(9, 5))
ax.hist(finales, bins=bordes, density=True,
        alpha=0.6, label="Caminatas aleatorias")
ax.plot(x, p_teorica, linewidth=2, label="Solución analítica")
ax.set_xlim(-60, 60)
ax.set_xlabel("Posición (µm)")
ax.set_ylabel("Densidad (1/µm)")
ax.set_title("Comparación con el modelo continuo")
ax.legend()
plt.tight_layout()
plt.show()

La distribución exacta de la caminata discreta es binomial. La distribución gaussiana es una aproximación que adquiere validez bajo las condiciones del límite difusivo. Las diferencias pueden proceder tanto de la naturaleza discreta del modelo como de las fluctuaciones estadísticas de la muestra.

9.3. Comprobación del desplazamiento cuadrático medio

En la sección 7 demostramos que, para una caminata aleatoria unidimensional, el desplazamiento cuadrático medio satisface

Como todas las partículas parten del origen, podemos estimarlo mediante \(\operatorname{DCM}_M(t)=M^{-1}\sum_{k=1}^{M}X_k^2(t)\).

dcm = np.mean(posiciones**2, axis=0)
dcm_teorico = 2 * D * tiempos
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(tiempos, dcm, label="Experimento numérico")
ax.plot(tiempos, dcm_teorico, "--", label="Predicción teórica")
ax.set_xlabel("Tiempo (s)")
ax.set_ylabel("DCM (µm²)")
ax.set_title("Desplazamiento cuadrático medio")
ax.legend()
ax.grid(alpha=0.25)
plt.tight_layout()
plt.show()

La predicción teórica corresponde a una recta que pasa por el origen y cuya pendiente es \(2D\). El resultado experimental fluctuará alrededor de la predicción debido al número finito de partículas.

Estimación del coeficiente de difusión

Si ajustamos una recta \(\operatorname{DCM}_M(t)\simeq at+b\), la pendiente proporciona \(D_{\mathrm{estimado}}=a/2\).

pendiente, intercepto = np.polyfit(tiempos, dcm, 1)
D_estimado = pendiente / 2
error_relativo = abs(D_estimado - D) / D
print(f"D teórico: {D:.3f} µm²/s")
print(f"D estimado: {D_estimado:.3f} µm²/s")
print(f"Diferencia relativa: {100*error_relativo:.2f}%")

La estimación dependerá de la realización numérica. Los valores del desplazamiento cuadrático medio en instantes diferentes no son estadísticamente independientes, pues proceden de las mismas trayectorias. Por ello, un ajuste lineal ordinario proporciona una estimación útil, pero no una incertidumbre rigurosa.

El efecto del número de partículas

Podemos repetir el experimento con 100, 1 000 y 10 000 partículas. Al aumentar el tamaño de la muestra, suelen disminuir las fluctuaciones relativas de los promedios calculados. Sin embargo, aumentar el número de partículas no elimina las aproximaciones físicas del modelo.

9.4. Resolución numérica de la ecuación de difusión

Retomaremos el recipiente de longitud finita, sección transversal constante y paredes impermeables que estudiamos en las secciones 4 y 8. La concentración satisface

Aquí no representaremos partículas individuales: calcularemos directamente la concentración media de diferentes regiones del recipiente.

El método de volúmenes finitos

Dividiremos el recipiente en \(N_x\) celdas de anchura \(\Delta x=L/N_x\). Representaremos mediante \(c_i(t)\) la concentración media en la celda \(i\). El balance de materia exige

Al combinar ambos resultados, para las celdas interiores obtenemos

Aplicaremos Euler explícito con un paso temporal \(\Delta t\):

Para este esquema unidimensional, con \(D\) constante, la condición habitual de estabilidad es \(0<r\leq 1/2\). En las celdas extremas, los flujos externos serán nulos.

Un experimento numérico reproducible

Utilizaremos los parámetros de nuestro programa anterior: \(L=100\,\mu\mathrm m\), \(D=100\,\mu\mathrm m^2/\mathrm s\), \(N_x=201\), \(r=0.4\), \(T=30\,\mathrm s\), concentración de fondo \(c_b=0.1\,\mathrm{mol/m^3}\), amplitud \(c_a=0.9\,\mathrm{mol/m^3}\) y anchura inicial \(\sigma_0=5\,\mu\mathrm m\).

El programa calculará la evolución de este perfil y registrará distribuciones correspondientes a distintos instantes.

La conservación de la materia como prueba numérica

Cada flujo que sale de una celda entra en la vecina; al sumar los balances, los flujos interiores se cancelan. Con paredes impermeables la cantidad total debe permanecer constante:

Podemos verificarlo mediante el error relativo \(\varepsilon=|n_{\mathrm{final}}-n_{\mathrm{inicial}}|/n_{\mathrm{inicial}}\). Una implementación conservativa correcta debe producir un error muy pequeño, limitado fundamentalmente por los errores de redondeo.

Comparación con la solución analítica

Para las paredes impermeables obtuvimos en la sección 8 una solución en serie de cosenos. Podemos comprobar nuestro método con una condición inicial que contenga un solo modo: \(c(x,0)=\overline c+a_1\cos(\pi x/L)\). En tal caso,

Esta expresión ofrece una referencia analítica con la cual estudiar cómo cambia el error al modificar la cantidad de celdas o el paso temporal. Nos permitirá comprobar tanto la conservación como la convergencia del procedimiento.

9.5. Comparación de los tres métodos: ¿por qué podemos confiar en las simulaciones?

A lo largo de este laboratorio hemos desarrollado tres procedimientos para estudiar la difusión: caminatas aleatorias, solución analítica y volúmenes finitos. Para comparar sus predicciones, los aplicaremos al mismo problema físico y examinaremos sus semejanzas, diferencias y fuentes de error.

Un problema común

Consideremos un recipiente de longitud \(L=100\,\mu\mathrm m\), con paredes impermeables y un coeficiente constante \(D=100\,\mu\mathrm m^2/\mathrm s\). La concentración inicial es

Utilizaremos unidades relativas de concentración. Como demostramos en la sección 8, la solución analítica es

Dividimos el recipiente en 100 celdas de \(1\,\mu\mathrm m\) y simulamos 150 000 partículas independientes. Cada partícula da pasos hacia una celda vecina cada \(0.005\,\mathrm s\); cuando intenta atravesar una pared, permanece en la celda extrema. Así implementamos fronteras reflectantes compatibles con el flujo nulo. El método de volúmenes finitos utiliza la misma discretización y un parámetro \(r=D\Delta t/(\Delta x)^2=1/2\). Comparamos los resultados después de un segundo.

Comparación de los resultados

Comparación cuantitativa de solución analítica, volúmenes finitos y caminatas aleatorias
Figura 8. Comparación de los tres métodos para las mismas condiciones iniciales y de frontera, en \(t=1\,\mathrm s\). La solución analítica se promedia sobre cada celda para compararla coherentemente con los métodos discretos.

La solución analítica y la de volúmenes finitos son casi indistinguibles a la escala de la figura. La distribución obtenida con las partículas fluctúa alrededor de ambas curvas porque procede de una muestra finita. Los tres métodos describen el transporte de una región más concentrada hacia otra menos concentrada. En esta ejecución, el error máximo entre volúmenes finitos y la solución analítica promediada es aproximadamente \(8.8\times10^{-6}\) unidades relativas.

La coincidencia de los dos métodos computacionales tiene una explicación adicional: con \(r=1/2\), la evolución esperada de las partículas en esta red coincide con el esquema discreto de volúmenes finitos, incluidas las fronteras. Por eso, comparar ambos no equivale a realizar dos verificaciones independientes. La solución analítica continua constituye una referencia distinta.

Semejanzas y diferencias

MétodoVentajaLimitación
Caminatas aleatoriasPermiten examinar trayectorias y fluctuacionesLos promedios presentan incertidumbre por muestreo
Solución analíticaProporciona predicciones exactas dentro del modelo continuoNo siempre se puede obtener para problemas complejos
Volúmenes finitosSe adapta a condiciones y geometrías más generalesIntroduce errores de discretización y exige comprobar estabilidad y convergencia

¿Cómo comprobamos la convergencia?

Una simulación no adquiere credibilidad simplemente porque produzca una curva razonable. También debemos averiguar si la solución se aproxima a la del modelo matemático conforme refinamos la discretización.

Para estudiarlo resolvemos el mismo problema con 25, 50, 100 y 200 celdas. Elegimos el paso temporal de modo que \(r\leq0.4\) y terminamos todas las ejecuciones en \(t=1\,\mathrm s\). Comparamos la concentración calculada en cada celda con el promedio de la solución analítica sobre esa celda; así evitamos confundir el error numérico con una diferencia entre valores puntuales y promedios espaciales.

Convergencia del método de volúmenes finitos
Figura 9. Error máximo en función de la anchura de celda. La referencia discontinua es proporcional a \((\Delta x)^2\). Los datos muestran una disminución compatible con el orden esperado.
Celdas\(\Delta x\) (µm)Error máximo
254\(9.49\times10^{-5}\)
502\(2.44\times10^{-5}\)
1001\(6.18\times10^{-6}\)
2000.5\(1.54\times10^{-6}\)

Al reducir a la mitad la anchura de las celdas, el error disminuye aproximadamente por un factor de cuatro. Este resultado es coherente con la precisión de segundo orden espacial de nuestro esquema cuando el paso temporal disminuye proporcionalmente a \((\Delta x)^2\). Un estudio de convergencia más amplio permitiría examinar otras condiciones y distinguir con mayor detalle los errores espaciales y temporales.

De la verificación a la confianza

La verificación examina si el programa resuelve correctamente las ecuaciones elegidas. Incluye comprobar las leyes de conservación, contrastar soluciones conocidas y estudiar la convergencia. En una simulación estocástica también debemos cuantificar el efecto de utilizar un número finito de partículas.

La validación, en cambio, investiga si el modelo y sus hipótesis representan adecuadamente el fenómeno real. Una simulación puede resolver una ecuación con enorme precisión y, sin embargo, no describir bien un sistema biológico si sus hipótesis físicas son inadecuadas.

A medida que estudiemos sistemas más complejos —como membranas con conductancias iónicas, ecuaciones diferenciales acopladas y redes neuronales—, las soluciones analíticas serán menos frecuentes. Los métodos numéricos nos permitirán investigar esos sistemas, siempre que conozcamos sus aproximaciones y sometamos sus resultados a comprobaciones apropiadas.

Las simulaciones no sustituyen la comprensión física: se apoyan en ella, permiten comprobar sus consecuencias y amplían el conjunto de problemas que podemos estudiar.

9.6. Tres descripciones de un mismo fenómeno

Hemos utilizado diferentes procedimientos para estudiar la difusión. Aunque todos se refieren al transporte de materia, cada uno trabaja con variables y aproximaciones diferentes.

DescripciónVariable fundamentalCaracterística
Caminatas aleatoriasPosición de las partículasEstocástica
Solución analíticaConcentración o densidad de probabilidadContinua y determinista
Volúmenes finitosConcentración media de las celdasNumérica y determinista

La simulación de caminatas aleatorias permite estudiar las trayectorias individuales y calcular propiedades estadísticas del conjunto. La solución analítica describe la evolución continua bajo condiciones iniciales y de frontera específicas. El método de volúmenes finitos aproxima la solución cuando no disponemos de una expresión analítica sencilla o queremos estudiar condiciones más generales.

Nuestros dos experimentos principales, presentados en los apartados 9.1–9.4, no representan exactamente el mismo problema: uno parte de partículas localizadas en el origen y sin fronteras; el otro utiliza una distribución inicial de anchura finita en un recipiente cerrado. La subsección 9.5 resolvió esta dificultad al plantear expresamente un problema común a los tres procedimientos y comparar sus resultados bajo condiciones compatibles.

El alcance de nuestros resultados

Las simulaciones permiten explorar las consecuencias de nuestras hipótesis matemáticas, examinar sus predicciones y comprobar la implementación de los procedimientos numéricos. Sin embargo, la coincidencia entre una simulación y una solución analítica no demuestra, por sí sola, que el modelo represente correctamente un sistema biológico real. Para establecer esa correspondencia necesitamos contrastar sus predicciones con observaciones experimentales y determinar si las hipótesis adoptadas son apropiadas.

Esta distinción entre modelo, resolución matemática y evidencia experimental será fundamental a lo largo de nuestro curso.

10. De la difusión ideal a la membrana neuronal

Hemos recorrido un camino que comenzó con una pregunta aparentemente sencilla: ¿cómo puede el movimiento desordenado de las partículas producir un transporte neto de materia? Para responderla, introdujimos magnitudes macroscópicas, dedujimos una ecuación de evolución y construimos una descripción estadística basada en caminatas aleatorias. Después estudiamos soluciones analíticas y utilizamos simulaciones para contrastar los distintos procedimientos.

Sin embargo, el propósito de este laboratorio no era estudiar la difusión como un fenómeno aislado. Nuestro interés está en comprender cómo se transportan las sustancias en el sistema nervioso y, particularmente, qué ocurre cuando ese transporte involucra una membrana neuronal.

Ha llegado el momento de examinar qué nos permiten explicar nuestros modelos, qué aspectos dejan fuera y cuáles serán las nuevas preguntas físicas que tendremos que abordar.

10.1. Lo que hemos aprendido de la difusión

La primera ley de Fick establece una relación entre el flujo difusivo y las variaciones espaciales de la concentración. En el modelo unidimensional que hemos considerado, esta relación se expresa como

\[J=-D\frac{\partial c}{\partial x}.\]

El signo negativo indica que, en las condiciones de validez de esta ley, el transporte neto ocurre en sentido contrario al gradiente de concentración. La ecuación no afirma que cada partícula se mueva hacia las regiones menos concentradas: las partículas continúan desplazándose de manera irregular y pueden atravesar una misma superficie en ambos sentidos. El flujo representa el resultado estadístico de esos movimientos.

Al combinar la primera ley de Fick con el principio de conservación de la materia, y suponiendo un coeficiente de difusión constante, obtuvimos la ecuación de difusión:

\[\frac{\partial c}{\partial t}=D\frac{\partial^2c}{\partial x^2}.\]

Esta ecuación permite predecir cómo evoluciona una distribución de concentraciones cuando conocemos su estado inicial y las condiciones que se cumplen en las fronteras. También explica por qué las variaciones espaciales tienden a suavizarse con el tiempo.

La descripción estadística aportó otra perspectiva. A partir de una caminata aleatoria simétrica, cuyos pasos son independientes y carecen de una dirección privilegiada, encontramos una relación entre el coeficiente de difusión, la longitud de los pasos y el intervalo temporal. En una dimensión, el desplazamiento cuadrático medio satisface

\[\left\langle [X(t)-X(0)]^2\right\rangle=2Dt.\]

Por tanto, la escala espacial característica de la dispersión aumenta como la raíz cuadrada del tiempo, no de manera lineal. Esta propiedad resulta importante cuando comparamos procesos de transporte que actúan a distancias distintas.

Finalmente, las soluciones analíticas y los procedimientos numéricos nos permitieron estudiar un mismo fenómeno mediante representaciones complementarias. Su comparación mostró cómo podemos verificar una simulación en problemas cuya solución conocemos, sin confundir esa verificación matemática con la validación experimental de un modelo físico.

10.2. Los límites de nuestro modelo

Todas las conclusiones anteriores se obtuvieron bajo determinadas hipótesis. Supusimos un medio homogéneo, un coeficiente de difusión constante, ausencia de transporte macroscópico del fluido y una sustancia que no se produce ni se consume mediante reacciones químicas. En el modelo de caminata aleatoria, además, consideramos pasos independientes y probabilidades de desplazamiento simétricas.

Estas aproximaciones nos permitieron identificar los mecanismos fundamentales sin tener que describir todos los detalles microscópicos. No obstante, debemos evitar interpretar la ecuación obtenida como una ley que describa cualquier forma de transporte en cualquier medio.

En un medio heterogéneo, por ejemplo, el coeficiente de difusión puede depender de la posición. En ese caso, no podemos extraerlo de la derivada espacial y la ecuación adopta, en una dimensión, la forma

\[\frac{\partial c}{\partial t}=\frac{\partial}{\partial x}\!\left(D(x)\frac{\partial c}{\partial x}\right),\]

si mantenemos la ley de Fick como relación constitutiva y seguimos suponiendo que no existen fuentes ni sumideros. Esta expresión muestra que modificar una hipótesis del modelo puede cambiar la ecuación que debemos resolver.

Si el fluido se mueve, también habrá que considerar el transporte de sustancia asociado a ese movimiento. Si ocurren reacciones químicas, el balance de materia deberá incluir los términos de producción y consumo correspondientes. Si las partículas interactúan de forma importante o su movilidad depende de la concentración, puede resultar necesario emplear relaciones de transporte más generales.

Las condiciones de frontera son igualmente decisivas. Un recipiente cerrado y una región que intercambia materia con su entorno pueden obedecer la misma ecuación de difusión en su interior y, aun así, presentar evoluciones muy diferentes.

Estas observaciones nos conducen a una lección metodológica: las ecuaciones no se aplican independientemente de las hipótesis utilizadas para construirlas. Antes de emplear un modelo en un sistema biológico, debemos examinar cuáles de esas hipótesis siguen siendo razonables.

10.3. La membrana como frontera selectiva

Una neurona no es simplemente una región de líquido en contacto directo con otra. Su interior y el medio extracelular están separados por una membrana plasmática, constituida principalmente por una bicapa de lípidos en la que se encuentran incorporadas distintas proteínas.

La presencia de esta estructura modifica profundamente el problema del transporte. Algunas sustancias pueden atravesar la bicapa lipídica con relativa facilidad, mientras que otras encuentran una barrera importante. Los iones, debido a su carga eléctrica y a las interacciones que mantienen con el agua y con su entorno, generalmente requieren vías específicas para atravesarla en cantidades fisiológicamente significativas.

Entre esas vías se encuentran los canales iónicos, que proporcionan recorridos selectivos para determinadas especies, y los transportadores, que permiten otros mecanismos de intercambio. Algunas proteínas transportadoras aprovechan fuentes de energía para mantener diferencias de concentración que no podrían sostenerse mediante transporte pasivo exclusivamente.

Por tanto, no basta con conocer la concentración de una sustancia a ambos lados de la membrana. También necesitamos saber si existe una vía que le permita cruzarla, cuáles son las propiedades de esa vía y en qué condiciones se encuentra disponible.

En un modelo macroscópico sencillo, esta dificultad puede representarse mediante una condición de frontera que relacione el flujo a través de la membrana con las variables físicas a ambos lados de ella. Sin embargo, esa relación debe elegirse de acuerdo con el mecanismo considerado: no será necesariamente la misma para una molécula que atraviesa la bicapa, un ion que pasa por un canal o una sustancia transportada mediante un proceso que consume energía.

La selectividad de la membrana introduce, así, una diferencia esencial respecto de nuestro recipiente idealizado. El transporte ya no depende únicamente de cómo se distribuye la sustancia en el espacio, sino también de las propiedades de la frontera que separa los dos medios.

10.4. Cuando las partículas tienen carga eléctrica

Hasta ahora hemos descrito partículas cuyo transporte depende solamente de las variaciones de concentración. Pero muchas de las sustancias más importantes para la actividad neuronal son iones, es decir, partículas que poseen carga eléctrica. Entre ellas se encuentran los iones de sodio, potasio, calcio y cloruro.

Un ion situado cerca de una membrana puede experimentar simultáneamente dos tendencias. Por un lado, su movimiento térmico da lugar a un transporte difusivo asociado a las diferencias de concentración. Por otro, si existe un campo eléctrico, su carga hace que experimente una fuerza eléctrica.

La fuerza sobre una partícula de carga \(q\), situada en un campo eléctrico \(\mathbf E\), está dada por

\[\mathbf F_{\mathrm{el}}=q\mathbf E.\]

El sentido de esta fuerza depende del signo de la carga. Los iones positivos y negativos responden en direcciones opuestas ante un mismo campo eléctrico. Por eso, una diferencia de potencial entre ambos lados de la membrana puede favorecer el desplazamiento de ciertos iones y oponerse al de otros.

Es importante distinguir las fuerzas que actúan sobre una partícula del flujo neto que aparece al considerar muchas partículas. El flujo total depende del comportamiento colectivo y de las propiedades del medio y de la membrana. No podemos obtenerlo simplemente sumando una concentración y una fuerza eléctrica, pues se trata de magnitudes físicas distintas.

En determinadas condiciones, la tendencia difusiva y la contribución del campo eléctrico pueden compensarse. Entonces es posible que el flujo neto de una especie iónica sea nulo, aunque sus partículas continúen en movimiento. Este estado se relaciona con el concepto de equilibrio electroquímico.

La relación cuantitativa entre la concentración, la carga eléctrica y el potencial se desarrollará cuando estudiemos los potenciales de equilibrio. Más adelante incorporaremos también las ecuaciones de transporte que describen conjuntamente la difusión y los efectos eléctricos. Por ahora, lo fundamental es reconocer por qué nuestro modelo puramente difusivo resulta insuficiente para describir el comportamiento de los iones a través de una membrana neuronal.

10.5. De los modelos sencillos a los sistemas complejos

El análisis de la difusión nos ha proporcionado un ejemplo de cómo se construye un modelo físico: identificamos las variables pertinentes, formulamos hipótesis, establecemos una relación de transporte y aplicamos un principio de conservación. Después resolvimos las ecuaciones en algunos casos y comprobamos mediante simulaciones que los procedimientos numéricos reproducen las predicciones conocidas.

Este modo de trabajo seguirá siendo útil cuando estudiemos la actividad eléctrica neuronal. La membrana no solo separa medios con distintas composiciones químicas; también puede mantener una diferencia de potencial eléctrico y modificar su estado cuando cambian las corrientes que la atraviesan.

Para comprender estos fenómenos tendremos que combinar los mecanismos de transporte iónico con las propiedades eléctricas de la membrana. Los circuitos equivalentes nos permitirán representar algunas de esas relaciones mediante elementos como capacitancias y conductancias. Las ecuaciones diferenciales describirán su evolución temporal y las simulaciones numéricas nos ayudarán a investigar el comportamiento de los sistemas cuando las variables estén acopladas.

Conforme aumente la complejidad del modelo, será menos frecuente encontrar soluciones analíticas completas. Esto no significa que debamos renunciar a la comprensión: podremos estudiar casos límite, comprobar dimensiones y leyes de conservación, analizar cómo influyen los parámetros y comparar los resultados numéricos con observaciones experimentales.

También deberemos mantener una distinción que ha acompañado todo nuestro laboratorio. Una simulación puede ayudarnos a comprender las consecuencias de un conjunto de ecuaciones; determinar si esas ecuaciones describen adecuadamente una neurona exige contrastarlas con la evidencia biológica.

De esta manera, la difusión constituye nuestro primer ejemplo de una estrategia que utilizaremos reiteradamente: comenzar con un fenómeno sencillo, construir un modelo que permita explicar sus mecanismos, comprobar sus predicciones y ampliar gradualmente la descripción cuando aparezcan nuevas interacciones.

Hemos comprendido cómo un movimiento microscópico irregular puede producir leyes macroscópicas de transporte. El siguiente problema será incorporar las propiedades eléctricas de las partículas y de la membrana que deben atravesar.

Si los iones tienden a difundirse siguiendo sus gradientes de concentración, pero al mismo tiempo experimentan fuerzas eléctricas, ¿cómo puede establecerse un equilibrio entre ambas tendencias y qué consecuencias tiene para el potencial eléctrico de la membrana neuronal?

Bibliografía del laboratorio

Las siguientes obras permiten profundizar en los fundamentos físicos, estadísticos y numéricos de los modelos desarrollados en este laboratorio.

Crank, J. (1975). The Mathematics of Diffusion, 2.ª edición. Clarendon Press, Oxford.

Berg, H. C. (1993). Random Walks in Biology, edición ampliada. Princeton University Press.

LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems. SIAM. DOI: 10.1137/1.9780898717839.