Fenômenos estocásticos em migração celular: teoria, experimentos e simulações

Stochastic phenomena in single cell migration: theory, experiments, and simulations

Resumos

Migração celular é um ingrediente importante no desenvolvimento de seres multicelulares e em processos como resposta imunológica ou metástases em câncer. Apresentamos uma revisão sobre modelagem de migração celular de célula isolada sobre substratos planos. Reanálise de dados experimentais mostraram um comportamento que se desvia do modelo canônico, o modelo de Langevin para uma partícula com um movimento Browniano. A proposição de um modelo semiempírico que prevê um comportamento difusivo em escalas de tempo muito curtas, ajusta os dados experimentais mas coloca em cheque a definição usual de velocidade e o correspondente protocolo de medidas. A solução é apresentada sob a forma de um modelo anisotrópico para a migração celular com uma variável interna representando a quebra de simetria espacial da célula em migração. A previsão teórica possibilita a definição de unidades naturais do sistema que resultam no colapso de curvas teóricas e experimentais em uma família de curvas com um só parâmetro. Adicionalmente, discutimos um modelo de simulações para células em três dimensões, que foi validado por comparações quantitativas com dados experimentais, possibilitadas pelo modelo teórico proposto. Estes resultados podem agora ser aplicados para situações mais complexas, que implicam em interações não lineares e requerem soluções numéricas.

Palavras-chave
Migração celular; Equação de Fürth modificada; Processos de Ornstein-Uhlenbeck; CompuCel3D

Cell migration plays important roles in the development of multicellular organisms and processes as immune response and in cancer metastasis. We present here a review on cell migration modelling of single cells on a flat substrate. Re-analyses of experimental data evinced a behavior that deviates from the canonical Langevin model of a particle immersed on a viscous fluid that presents a Brownian persistent motion. The proposition of a semi-empirical model yields a diffusive behavior in short-time intervals and fits the experimental data, but challenges velocity definition and the corresponding measurement protocol. The solution is presented on the form of an anisotropic model for cell migration, that considers an additional, internal variable, that accounts for the spatial symmetry break of a migrating cell. The theoretical results yield the proposition of natural units for the problem that allows the collapse of theoretical and experimental curves onto a single parameter family of curves. We also discuss a simulation model for three dimensional cells migrating on a flat substrate. The simulation results are quantitatively validated using the theoretical results. The simulation model is now ready to be used in investigating cell migration in more complex environments, that imply non-linear interactions and require numerical solutions.

Keywords
Cell migration; Modified Fürth equation; Ornstein-Uhlenbeck processes; CompuCel3D

1. Introdução

Complexidade é difícil de definir. Há muita confusão e talvez seja assim porque as pessoas se referem como complexos a sistemas muito diferentes. Há diferenças importantes entre uma visão experimental e uma visão teórica do que significa a complexidade. Para eu explicar essa diferença, preciso antes dar um passo atrás e discutir um pouco qual é, na minha concepção, uma definição operacional para Física.

Suponho ser sem controvérsias que Física utiliza matemática para descrever fenômenos. Mas a Natureza, o palco dos fenômenos, não é composta de números. Assim, há que traduzir os fenômenos em números e isto se faz através de medidas. Em outras palavras, protocolos de medidas são os códigos que traduzem os valores assumidos pelas variáveis de um sistema matemático formal para as características de um sistema natural. Por exemplo, quando usamos as equações de Newton para descrever o lançamento de um projétil, associamos os resultados do modelo às medidas de distância, altura, tempo e velocidade.

Medidas, por outro lado, nem sempre são simples e diretas. Há sutilezas importantes. Excetuando talvez as grandezas cujas medidas são as fundamentais, aquelas associadas às unidades fundamentais do sistema de unidades, todas as outras medidas dependem de um modelo teórico. No caso do Sistema Internacional (SI), as grandezas fundamentais são comprimento, tempo, massa, corrente elétrica, temperatura termodinâmica, quantidade de matéria, e intensidade luminosa, cujas unidades são, respectivamente, metro, segundo, quilograma, Ampère, Kelvin, mol e candela.

A velocidade de uma bola, por exemplo, não é medida por uma unidade fundamental, mas por uma composição delas, m/s. A velocidade média de um móvel em um determinado intervalo de tempo Δt, é definida como a razão entre o deslocamento $\mathrm{\Delta }\overrightarrow{r}$ e o intervalo de tempo. A velocidade instantânea, por outro lado, é definida pelo limite desta razão quando o intervalo de tempo vai a zero. Se decidirmos medir velocidade instantânea usando a medida de deslocamentos e intervalos de tempo, uma medida robusta da velocidade requer que a razão $\mathrm{\Delta }\overrightarrow{r}/\mathrm{\Delta }t$ convirja para um valor finito se Δt se aproxima de zero. Em termos práticos, precisamos garantir que haja um limite de tamanho para o intervalo de tempo tal que o valor dessa razão não mudará para intervalos menores. Se existir tal limite e ele for acessível aos instrumentos de medida que utilizamos, teremos um protocolo robusto para a medida de velocidade.

Retornando à ideia de definição de complexidade, podemos definir o grau de complexidade pelo número de variáveis que precisamos conhecer para completamente prever a evolução do sistema. Então, quais sistemas poderiam ser classificado de um sistema complexo? A minha resposta é que depende do que se define como sistema. Em primeiro lugar é preciso diferenciar se estamos falando de sistemas naturais ou sistemas matemáticos formais.

No primeiro caso, arrisco a dizer que só existe um sistema natural, que é o Universo. Dado tempo suficientemente grande, uma dada parte do universo pode ter seu estado alterado pelo que acontece em outro ponto. Um exemplo literalmente bombástico é a colisão do meteoro com a Terra há 65 milhões de anos que provocou a extinção da maioria das espécies de dinossauros e possibilitou que a evolução dos mamíferos fosse tal que nós pudéssemos surgir como espécie: a criação e trajetória do meteoro em lugares distantes da nossa galáxia determinou a trajetória seguida pelo Ecossistema terrestre. Assim, para descrevermos o Ecossistema terrestre durante todo o período de sua existência, teríamos de considerar a dinâmica da formação de meteoros em regiões distantes na nossa galáxia. O que significa um número muito grande de variáveis levando à conclusão que o sistema natural representado pelo Universo é complexo. Observe que, em sistemas naturais, inevitavelmente há não linearidades se considerarmos sistemas muito extensos e tempos muito longos.

Por outro lado, se nossa ambição não for tão grande e limitamos o tempo de validade de nossa descrição de uma parte também limitada do Universo, o número de variáveis que precisamos considerar diminui consideravelmente. Mas, neste caso, estaremos fazendo aproximações, o que necessariamente implica em hipóteses e modelagem. Em Física, modelagem passa necessariamente pela utilização de técnicas matemáticas e protocolos de medidas. Considerando então sistemas matemáticos formais podemos definir que sua complexidade é de alguma forma medida pelo número de variáveis que precisamos conhecer para prevermos o que acontece. Assim, o modelo de gás ideal termodinâmico, descrito em termos de 3 variáveis – pressão, volume e temperatura – é um sistema simples. Já o modelo de gás ideal da mecânica estatística clássica é complexo, com 6×6,02×10²³ variáveis para um mol de gás, e pode apresentar turbulência. Observe que o fenômeno ao qual se aplicam os dois modelos – o gás monoatômico rarefeito a temperaturas ambientes – é o mesmo; a qualificação sobre a complexidade se aplica aos modelos, não ao fenômeno. Fenômenos acabam sempre por serem complexos.

A Física é um campo pródigo em modelos relativamente simples para sistemas naturais. Na verdade, excetuando-se raros casos (talvez a turbulência proposta por Kolmogorov se encaixe nesses casos), os sistemas que apresentam comportamentos que podem ser descritos por estes modelos matemáticos lineares e com um número limitado de variáveis é que puderam ser analiticamente resolvidos até meados do século XX. Como computadores não eram bem desenvolvidos e disseminados como hoje, técnicas numéricas também não podiam ser utilizadas em todo o seu potencial. Assim, a mecânica estatística focava em casos de equilíbrio termodinâmico, mecânica de fluídos em fluxos lineares e assim por diante.

Fenômenos biológicos, por outro lado, não são em geral passíveis de uma modelagem matemática linear e com poucas variáveis. Na realidade, protocolos robustos para medidas em sistemas biológicos são especialmente difíceis: há frequentemente variáveis de confusão. Tipicamente, uma variável de confusão é associada a uma característica do sistema biológico que não foi levada em conta. Por exemplo, a fase do ciclo celular das células em um teste de resposta ao tratamento por uma droga: em diferentes fases as células podem responder com diferentes intensidades, aumentando a variância de uma medida.

O número de variáveis de confusão em sistemas celulares pode crescer com a diversidade genética do conjunto de células sob observação. Dependendo do que se está observando, isso equivale a dizer que o número de variáveis necessárias para a previsão de um fenômeno tende a crescer. Em outras palavras, o modelo usado, para ser eficiente, precisar aumentar o número de variáveis e, portanto, sua complexidade.

Resumindo, complexidade é um conceito útil quando aplicado ao modelo, mais do que ao fenômeno em si. Em fenômenos biológicos, a complexidade do modelo depende do objetivo (o modelo mais simples para descrever a vida de um mamífero é prever que todos o mamífero vivo vai morrer: absolutamente correto, mas não traz informação útil). Para construir um modelo para fenômenos biológicos é necessário a definição de medidas robustas. Modelos para fenômenos biológicos que trazem informação precisam lidar com variáveis de confusão, o que pode torná-los mais complexos.

Neste capítulo vamos tratar da modelagem de migração celular, um fenômeno biológico. Na seção seguinte, há uma discussão sobre a relevância de migração celular, e uma pequena revisão do modelo comumente aplicado à migração celular – a equação de Fürth – apresentando seu sucesso e limitações. Na seção 3 fazemos uma rápida revisão de biologia e, na seção 4, discutimos a consequência teórica e experimental da divergência dos dados experimentais em relação à equação de Fürth e discutimos um modelo teórico que resolve e propõe protocolos de medidas cinéticas suficientemente robustas para a comparação de resultados obtidos para diferentes células, com diferentes montagens experimentais e para a tradução para dados obtidos de simulações numéricas. Na seção 6 revisamos um modelo de simulação de migração celular, com células extensas e em 3 dimensões. Finalmente, na seção 6 concluimos e discutimos perspectivas e então voltamos à discussão sobre a complexidade em migração celular.

2. O Fenômeno

Em muitos processos em organismos pluricelulares, a migação celular desempenha um papel central. Podemos citar por exemplo o desenvolvimento do embrião [¹], em processos de recuperação de ferimentos [², ³] e na resposta inflamatória [⁴]. Por exemplo, uma modelagem eficiente da migração celular pode contribuir para a compreensão dos fatores que envolvem o recrutamento de fibroblastos para o sítio de inflamação nos pulmões, para o caso atual da COVID-19, ou para a engenharia de órgãos e tecidos celulares visando aplicações médicas e tecnológicas [⁵].

A migração de células eucariotas é um fenômeno para o qual podemos identificar diferentes escalas. Para uma revisão veja [⁶]. No interior da célula, as interações acontecem na escala molecular, embora caiba ressaltar que moléculas bioquímicas podem ser muito grandes, contendo milhares de átomos e podendo apresentar diferentes conformações espaciais. Há a escala das organelas como a mitocôndria e as fibras de actina e microtúbulos, e há a escala da célula, que consiste o corpo que migra. A migração celular é um fenômeno de matéria ativa, no sentido que a energia necessária para o movimento é fornecida pela maquinaria celular, assim como os pontos de apoio e fixação em substratos para trocar momento com o ambiente. Um complicador adicional é que as várias escalas não são separáveis: a forma da célula pode modificar a produção de proteínas da célula, o que pode alterar sua forma e a comunicação com o meio externo, dando lugar a feedback loops que atravessam diferentes escalas [⁷, ⁸, ⁹, ¹⁰, ¹¹]. Sistemas naturais nunca são simples.

De qualquer maneira, um modelo não pode lidar com informação infinita. Há que reduzir essa informação, com base (i) nos objetivos do modelo e (ii) algum conhecimento prévio do que deve influir no fenômeno específico que queremos modelar. No entanto, imersos na abundância de dados e observações típicas de sistemas naturais, os cientistas não têm certeza a priori de quais dados devem ser levados em conta e quais descartar. A validação do modelo se dá a posteriori, com (i) a verificação de que o modelo descreve corretamente o que já se sabia e (ii) com a previsão ou descoberta de fatos novos, confirmados por experimentos. Resta observar que, como houve descarte de informação, necessariamente o modelo tem limites. Sim, modelar é uma arte que de exata não tem nada.

O primeiro passo em qualquer tentativa de modelagem de um fenômeno é necessariamente a observação e, depois, a proposição de medidas. A observação de migração celular pode se dar em muitos contextos, podendo ser in vivo ou in vitro. Normalmente o controle é mais fácil em observações in vitro.O ensaio clássico de medidas de migração celular consiste em colocar em uma placa de petri células imersas em um meio de cultura. Colocando essa montagem em um microscópio invertido, é possível observar as células migrando. Registrando a posição da célula em diferentes tempos, pode-se caracterizar cineticamente o movimento. Vale ressaltar que as células são corpos extensos e, portanto, é necessário descrever como a posição da célula é medida. Em geral toma-se o ponto mais central da projeção da célula sobre o substrato, ou o centro do núcleo celular. É uma escolha que depende da montagem do experimento.

Para experimentos deste tipo, em 1920 Fürth modelou migração das células únicas por meio do deslocamento quadrático médio (MSD, do inglês mean squared displacement), $⟨{|\overrightarrow{r}(t+\mathrm{\Delta }t)-\overrightarrow{r}(t)|}^2⟩$, que dá a média do módulo quadrático do deslocamento apresentado pelas células como função do intervalo de tempo Δt usado para a medida do deslocamento [¹²]. A média é tomada sobre todos os tempos t e sobre diferentes células. A equação de Fürth para o MSD é

(1) $$MSD=⟨{|\overrightarrow{r}(t+\mathrm{\Delta }t)-\overrightarrow{r}(t)|}^2⟩=4D(\mathrm{\Delta }t-P(1-e^{-\mathrm{\Delta }t/P})),$$

onde P é conhecido como o tempo de persistência e D, o coeficiente de difusão. Observe que para pequenos Δt, $MSD\sim \frac{2D}{P}\mathrm{\Delta }{t}^2$, reproduzindo um comportamento balístico em que o deslocamento é proporcional a Δt, enquanto que para grandes intervalos de tempo, MSD∼4DΔt, mostrando um comportamento difusivo. A escala de tempo para a qual o comportamento balístico persiste é dada por P, que marca a transição entre os dois regimes.

Pela equação de Fürth, as células descreveriam um movimento balístico para intervalos de tempos pequenos e difusivo com coeficiente D em intervalos maiores que o tempo de persistência P, parâmetros típicos de cada célula. Reforçando a modelagem, a Equação (1) pode ser obtida como a solução de um problema elegantemente colocado por Langevin em 1908 [¹³] e resolve um problema de longa data para o movimento Browniano [¹⁴].

O comportamento balístico para intervalos de tempo curtos resolve um problema conceitual importante do movimento Browniano para corpos pontuais: quando o deslocamento é proporcional a Δt e existe o limite ${lim}_{\mathrm{\Delta }t\to 0}\frac{\mathrm{\Delta }\overrightarrow{r}}{\mathrm{\Delta }t}$, velocidade é um conceito bem definido e é possível a definição de protocolos reprodutíveis de medidas.

No entanto, para o experimento de Brown de 1827, com partículas liberadas de grão de pólen em água, tal intervalo parecia não existir: à medida que Δt diminui, o módulo do deslocamento do centro de massa da partícula é proporcional a $\sqrt{\mathrm{\Delta }t}$, de tal maneira que ${lim}_{\mathrm{\Delta }t\to 0}\mathrm{\Delta }\overrightarrow{r}/\mathrm{\Delta }t$ não existe. A solução é então caracterizar o movimento pela média da razão ${|\mathrm{\Delta }\overrightarrow{r}|}^2/\mathrm{\Delta }t$, que define o coeficiente de difusão D deste movimento [¹⁵]. Na realidade, o movimento Browniano compartilha muitas características de fenômenos difusivos. Por outro lado, se velocidade instantânea não é uma quantidade bem definida, aceleração também não o é e a aplicação das leis de Newton entra em cheque.

Difusão é um fenômeno muito comum para que se possa conviver sem uma teoria robustamente embasada nas leis fundamentais da física. Ou pelo menos Einstein [¹⁵], Smoluchowski [¹⁶], Langevin [¹³] e Wiener [¹⁷] não puderam conviver com este problema. A solução é considerar uma partícula clássica, que segue as leis de Newton, imersa em um líquido. Se a partícula se move, o líquido exerce sobre ela uma força proporcional à velocidade, mas em sentido contrário, isto é, um atrito de Stokes, que dissipa energia. Ao mesmo tempo, as moléculas do líquido continuamente colidem com a partícula cedendo-lhe momento linear. A distribuição do impulso cedido nas colisões segue uma distribuição Gaussiana isotrópica, com média zero e um desvio padrão que depende do valor de velocidade quadrática média das moléculas, isto é, depende da temperatura do líquido. Uma colisão dessas é modelada como acontecendo instantaneamente, de tal maneira que induz mudanças descontínuas na velocidade. O ingrediente desta teoria que salva o conceito de velocidade instantânea é que há um tempo típico, pequeno, mas finito, entre duas colisões. Para intervalos de tempo menores que esse tempo entre colisões, a velocidade é bem definida (no sentido que que ${lim}_{\mathrm{\Delta }t\to 0}\mathrm{\Delta }\overrightarrow{r}/\mathrm{\Delta }t$ existe). Assim, o resultado é que a partícula recebe impulsos aleatoriamente distribuídos e perde velocidade no período entre colisões devido ao atrito de Stokes. O ganho de energia pelas colisões é compensado pela perda de energia devido ao atrito (teorema de flutuação-dissipação) e a partícula descreve uma trajetória estacionária, com energia cinética média igual à energia cinética média das moléculas do fluido. O modelo matemático correspondente a este problema foi proposto por Langevin [¹³] e pode ser escrito como

(2) $${\displaystyle \frac{d\overrightarrow{v}}{dt}}=-\gamma \overrightarrow{v}+\overrightarrow{\xi }(t)$$ $${\displaystyle \frac{d\overrightarrow{r}}{dt}}=\overrightarrow{v},$$

onde $\overrightarrow{r}$ e $\overrightarrow{v}$ são, respectivamente, a posição e a velocidade da partícula, γ representa a fração de velocidade perdida em dt devido à resistência do fluído e $\overrightarrow{\xi }(t)$ representa um ruído branco Gaussiano.

A solução estacionária deste problema leva a uma curva de deslocamento quadrático médio como função do intervalo de tempo que tem a mesma forma da Equação (1): apresenta dois comportamentos, tendo a transição entre eles de forma suave. Se Δt é bem menor que o tempo de colisões com as moléculas do líquido o movimento é balístico, isto é, $⟨{|\mathrm{\Delta }\overrightarrow{r}|}^2⟩\sim \mathrm{\Delta }{t}^2$ e se Δt é bem maior que este tempo entre colisões o movimento é difusivo, com $⟨{|\mathrm{\Delta }\overrightarrow{r}|}^2⟩\sim \mathrm{\Delta }t$. Assim, para tempos de observação maiores que o tempo de colisões com as moléculas observaremos um movimento de difusão, mas a partícula tem uma velocidade instantânea bem definida, que pode ser modelada usando mecânica Newtoniana. A saída para o problema de Brown é que no problema das partículas do pólen, o intervalo de tempo entre colisões não seria acessível aos aparelhos de medidas utilizados. Haveria ainda que entender as mudanças bruscas de velocidade devido às colisões e como devemos tratar matematicamente este problema. Para isso a contribuição de Wiener e mais tarde de Ito e Stratonovich foram essenciais para o desenvolvimento da área. A história de como a teoria sobre o movimento Browniano evolui é bem apresentada por Genthom [¹⁴], enquanto que uma revisão mais matemática e com generalizações interessantes é apresentada em [¹⁸].

Voltando às células, o experimento de Fürth de 1920 possibilitou a medida de intervalos de tempo pequenos e a explicação da Equação de Fürth foi tida como sendo a solução de um problema de Langevin. Só que as células não se comportam bem assim. Recentemente, Thomas e colaboradores [¹⁹] analisaram os dados de experimentos mais recentes de 5 laboratórios diferentes, considerando 12 montagens experimentais diversas (tipo celular, meio de cultura, substrato diferentes) e mostraram um desvio no comportamento do MSD previsto para pequenos intervalos de tempo. Os dados mostram uma fase difusiva adcional para intervalos de tempo muito pequenos, o que significa que migração celular não pode ser caracterizada por um termo de velocidade instantânea, como partículas de um fluido. A equação usada para fitar o MSD das células é a equação de Fürth modificada [¹⁹]:

(3) $$MS{D}_{F\ddot{u}rthModificada}=2D\left[\frac{\mathrm{\Delta }t}{(1-S)}-P(1-e^{-\mathrm{\Delta }t/P})\right],$$

onde 0≤S≤1 é um parâmetro adimensional. A equação de Fürth modificada apresenta 3 regimes: para intervalos de tempo muito curtos, Δt≪SP, o comportamento é difusivo, com $MS{D}_{F\ddot{u}rthModificada}\sim \frac{2DS}{1-S}\mathrm{\Delta }t$, se SP < Δt < P o deslocamento quadrático médio apresenta um comportamento tipo balístico, dependendo de (Δt)^α, com α > 1 e crescente, e se Δt≫P o comportamento volta ser balístico, com $MS{D}_{F\ddot{u}rthModificada}\sim \frac{2D}{1-S}\mathrm{\Delta }t$. Por outro lado, se S=0, voltamos a ter a equação original de Fürth.

A equação de Fürth modificada pode ser re-escalada pela definição de unidades de comprimento e tempo usando os parâmetros P, DeS da equação, de tal maneira que podemos definir tempo e distância em unidades naturais do sistema como

(4a) $$\tau \equiv \frac{t}{P}$$

e

(4b) $$\overrightarrow{\rho }\equiv \frac{\overrightarrow{r}}{\sqrt{2DP/(1-S)}},$$

de tal maneira que a Equação (3) fica

(5) $$⟨{|\mathrm{\Delta }\overrightarrow{\rho }|}^2⟩=\mathrm{\Delta }\tau -(1-S)(1-e^{-\mathrm{\Delta }\tau }).$$

Neste caso, resta apenas Scomo um parâmetro adimensional. Assim o ajuste das medidas de deslocamento quadrático médio pela equação (5) pode ser visto como um protocolo de medida das grandezas P, D e S, como detalhado na Ref. [¹⁹]. Observe que, neste caso, o deslocamento quadrático médio depende apenas do único parâmetro restante, S. A Figura 1 mostra a equação de Fürth modificada para diferentes valores de S.

Figura 1
Gráfico log-log da Equação de Fürth modificada, em unidades naturais, para diferentes valores de S. A Equação de Fürth original corresponde ao caso em que S=0, mostrado pela linha tracejada.

A Figura 2 mostra os resultados da re-análise de Thomas e colaboradores para os dados experimentais, com o ajuste usando a equação (5). Observe também que a curva para S=0, em pontilhado, não descreve bem os dados experimentais para Δτ pequenos.

Figura 2
MSD em unidades naturais para 9 experimentos de 4 laboratórios diferentes [²⁰, ²¹, ²², ²³], mostrando o colapso dos dados sobre a equação de Fürth modificada. Ajustes como detalhados na Ref. [¹⁹]. A linha tracejada corresponde à equação de Fürth original (equação (1)), quando S=0. Nos insets, o mesmo para os dados em unidades do laboratório.

Uma segunda medida experimental que desvia do comportamento previsto diz respeito à autocorrelação da velocidade. Para sistemas já no estado estacionária, a autocorrelação da velocidade pode ser calculada a partir da segunda derivada do MSD [²⁴]. A Figura 3 mostra a função de autocorrelação de velocidade como função do intervalo de tempo entre as duas medidas de velocidades. O cálculo desta função foi feito a partir das trajetórias das células, usando diferentes valores de intervao de tempo δ para estimar velocidade como a razão entre deslocamento e intervalo de tempo. Nesta figura, para valores de δ < S, a função de correlação decai rapidamente à medida que Δτ diminui, mostrando que para Δτ e δ pequenos a velocidade se comporta como se a trajetória fosse difusiva.

Figura 3
Função de autocorrelação de velocidades para os dados de Wu e colaboradores [²²], tomando estimativas de velocidades considerando intervalos de tempo finitos δ. Quando δ < S, a função de autocorrelação decai quando Δτ diminui, indicando um regime difusivo para intervalos de tempo curtos [¹⁹].

Se, por um lado, a proposição da equação de Fürth modificada resolve o problema empírico da descrição da migração celular para intervalos de tempo pequenos, por outro lado ela traz três problemas que não podem ser deixados de lado. O primeiro, relaciona-se ao fato de não se poder aplicar a mecânica Newtoniana no caso de não se poder definir velocidade instantânea, o ponto tão discutido no início do século XX. O segundo, diz respeito ao modelo cinético que pode dar lugar à equação de Fürth modificada. O terceiro tem a ver com o protocolo de medida para caracterizar a cinética das células.

Como verermos um pouco mais adiante, os três problemas são resolvidos primeiro por tratar a célula como um corpo extenso e deformável e não como uma partícula pontual e segundo, por considerar explicitamente que a célula em migração não apresenta simetria espacial, mas apresenta um eixo de polarização bem definido. Na seção seguinte, discutimos alguns aspectos sobre migração celular, necessários para justificar hipóteses para o modelo teórico que apresentamos na seção ⁴.

3. Um Pouco de Biologia

Como discutimos acima, para modelar apropriadamente um fenômeno precisamos definir o que significa “apropriadamente” explicitando nossos objetivos, e precisamos escolher quais aspectos vamos considerar e o que vamos descartar. Em ambos casos, é bom conhecer o sistema que estamos modelando. No que segue, vou discutir o que é relevante para o modelo estocástico apresentado na seção seguinte. Neste sentido, a coleção de informação abaixo já está filtrada, focando em um modelo em específico.

A migração de células mesenquimais é resultado de uma maquinaria finamente sintonizada e que recruta muitas das reações que acontecem no interior da célula. Assim sendo, a célula não é capaz de se dividir, ou diferenciar quando está migrando, pois sua organização está focada na migração. De fato, quando a célula começa o processo de síntese do DNA para uma eventual mitose, a morfologia da célula muda drasticamente, a célula adquire uma forma simétrica e interrompe a migração [²⁵].

Para que uma célula mesenquimal possa entrar em migração é preciso que haja gradientes de proteínas especializadas, as pequenas Rho-GTPases, das quais a proteína RAC e a proteína Rho são membros importantes. No fenótipo migratório há uma concentração de RAC no fronte da célula e de Rho na sua traseira. Estes gradientes ocorrem ativamente, isto é, as células gastam energia (ATP) para movimentar as proteínas, utilizando-se de transporte por meio dos microtúbulos, que já foram chamados de ‘as auto-estradas do transporte celular’ [²⁶]. Quanto mais fácil for para a célula o transporte destas proteínas, mais fácil será para a célula migrar, sendo interessante então que microtúbulos e as fibras de acto-miosina disponham-se no interior da célula preferencialmente no eixo frente-ré [²⁷, ²⁸].

O ambiente criado no fronte da célula, pelo acúmulo de proteínas adequadas, torna possível a dinâmica de ‘esteira’, típica da rede de fibra de actina. Acontece da maneira como ilustrada na Figura 4, que resume o modelo de Abercrombie [²⁹]. Os monômeros de actina polimerizam, formando fibras. Essa polimerização tem um sentido preferencial, porque os monômeros não são simétricos. Assim, em uma ponta os monômeros são acrescidos, mas na ponta oposta a actina é liberada, voltando a ser monômeros. A orientação das fibras assim formadas depende do ambiente local: devido ao acúmulo das proteínas adequadas há uma direção preferencial. No entanto, estas fibras podem se ramificar, sempre formando um ângulo de 74^∘, devido à forma das proteínas que possibilitam a ramificação e de como elas se conectam com a fibra original. Essa ramificação contínua acaba por gerar uma rede de fibras de actina que, ao crescer, empurra a membrana celular para o que vamos chamar de frente. O resultado disto é uma protrusão da membrana, muito fina em relação ao resto do corpo da célula. Essa parte da célula é denominada de lamelipódio. Concomitantemente com a protrusão, complexos proteicos transmembrânicos se criam e maturam, conectando-se à matriz extracelular e às fibras de actina, ancorando a rede de actina ao substrato. Esses pontos de adesão não se movem em relação às fibras ou ao substrato, mas a constante polimerização das fibras na frente da célula e a desmontagem das fibras nas pontas opostas faz com que, em relação ao centro de massa da célula, os pontos de adesão movam-se em direção da traseira. Na ré da célula o ambiente proteico é diferente e os pontos de adesão são desmontados. Além disso, uma rede de acto-miosina que une a frente da célula (que está sendo empurrada para frente devido à constante polimerização/despolimerização da rede de actina), juntamente com o córtex celular puxam a ré da célula para a frente. Com essa maquinaria a célula migra em uma dinâmica que pode ser comparada à dinâmica de esteira de um tanque de guerra. Para uma revisão, veja [⁶].

Figura 4
Dinâmica de esteira: migração celular sobre substratos planos. A frente da célula apresenta um lamelipódio, formado por uma rede densa de actina polimerizada. O eixo de polarização da célua, da ré à frente, garante o fluxo das proteínas necessárias para a ré e para a frente, criando gradientes que diferenciam o meio intracelular nessas localizações. A rede de actina prolonga-se para a frente e o ambiente celular propicia a formação de complexos moleculares transmembrânicos, formando os focos de adesão que ancoram a célula pelas adesões com o substrato e com o citoesqueleto. Simultaneamente, na ré da célula, o ambiente celular favorece a desmontagem das adesões e o córtex celular, juntamente com fibras de tensão, puxam a ré para frente. Os focos de adesão e os monômeros na rede de actina no lamelipódio não se movem em relação ao substrato, mas a contínua montagem à frente e desmontagem em regiões intermediárias do eixo ré-frente resulta em uma rede de actina cujo centro de massa avança para frente, como uma esteira de um tanque de guerra.

Resta ressaltar que a forma do lamelipódio em geral lembra um leque mais do que um braço, e o perímetro do lamelipódio é grande em relação ao perímetro da célula. Neste perímetro há ondulações rápidas em comparação com o tempo típico de persistência do lamelipódio. Essas ondulações somam-se, resultando em uma propulsão preferencialmente para a frente, mas há flutuações na direção perpendicular à polarização que têm um caráter puramente aleatório. Essas flutuações são causadas nas redes de actina, que por sua vez estão conectadas às fibras de acto-miosina e ao córtex celular. Como a medida de localização da célula envolve a localização do núcleo ou o centro geométrico da projeção da célula sobre o substrato, tal dinâmica têm como consequência flutuações na posição da célula, o que pode ser medida como uma difusão rápida, observável para intervalos de tempo muito menores que o tempo de persistência do lamelipódio.

A polimerização/despolimerização e a ramificação das fibras de actina têm um carácter estocástico. No entanto, quanto maior for a rede de actina em um dado local, maior a probabilidade de ela crescer neste local, resultando no que pode ser reconhecido como uma excitação local. Ao mesmo tempo, um crescimento da rede de actina diminui a disponibilidade de actina e das outras proteínas necessárias para a formação do lamelipódio em todo o resto da célula, e portanto, reduz a probabilidade de formação de outros frontes de migração globalmente. Observe que estocasticidade juntamente com excitação local e inibição global são as componentes necessárias para quebra espontânea de simetria [³⁰].

As consequências dessa complicada maquinaria celular relevantes para o modelo que vamos apresentar são as seguintes:

1. Formação do lamelipódio é estocástica e pode quebrar a simetria da célula espontaneamente;

2. O lamelipódio é instável, mas pode persistir por algum tempo;

3. A dinâmica de formação e reformação do lamelipódio, assim como o seu crescimento é resultado da maquinaria interna da célula, definindo a célula como um agente ativo;

4. Enquanto o lamelipódio persiste, a dinâmica de esteira define a velocidade da célula na direção da polarização;

5. A quebra de simetria implica na existência de um eixo de polarização da célula e o movimento da célula nesta direção é necessariamente diferente do movimento celular na direção perpendicular ao eixo de polarização.

6. Flutuações da membrana do lamelipódio acontecem simultaneamente em múltiplos pontos.

Na próxima seção discutimos um modelo que tem por objetivo resolver o problema trazido à luz pela re-análise apresentada por Thomas e colaboradores: se velocidade de migração celular não é uma quantidade bem definida e se o modelo de Langevin apoia-se em equações diferenciais para a velocidade, qual o modelo dinâmico que pode resultar na equação de Fürth modificada?

4. O Modelo Anisotrópico

O modelo de Langevin tem como resultado para o MSD a equação de Fürth, equação (1). A equação de Fürth modificada considera uma correção, dada pelo parâmetro adimensional S, que dá origem a um comportamento difusivo para intervalos de tempo muito curtos. Em princípio o que pode ter parecido uma pequena correção, coloca em cheque a formulação do modelo dinâmico, uma vez que um comportamento difusivo em intervalos de tempo curtos impossibilita a definição de velocidade e um protocolo robusto para a medida de velocidade.

O modelo que propomos [³¹] baseia-se em considerar um sistema anisotrópico, com uma direção de polarização como uma variável adicional do sistema. Observe que anisotropia é um aspecto importante em células em migração. Com ajuda da polarização, podemos supor duas dinâmicas diferentes: uma para a direção da polarização e a outra na direção perpendicular à polarização. Por outro lado, a polarização como um grau de liberdade adicional do sistema também precisa de uma equação de evolução.

Definimos como polarização o vetor $\overrightarrow{p}(t)=(cos\theta (t)$, sinθ(t)), onde θ(t) é um ângulo em relação ao referencial do laboratório, que segue a seguinte equação dinâmica:

(6) $$[\theta (t+\mathrm{\Delta }t)-\theta (t)]={\int }_t^{t+\mathrm{\Delta }t}{\beta }_{\perp }(t)\mathit{d}t,$$

onde β_⊥(t) é um ruído branco. Na direção da polarização, que modela a direção preferencial das fibras de actomiosina e microtúbulos, supomos que para um pequeno intervalo de tempo Δt a velocidade da célula obedece a

(7) $$v_{\parallel }{}^{final}(t)=\left[(1-\gamma \mathrm{\Delta }t)v_{\parallel }{}^{inicial}(t)+{\displaystyle {\int }_t^{t+\mathrm{\Delta }t}}{\xi }_{\parallel }(t)\mathit{d}t\right],$$

onde γ representa uma dissipação e ξ_∥(t) é também um ruído Gaussiano, com unidades adequadas, como veremos logo a seguir. $v_{\parallel }{}^{inicial}(t)$ and $v_{\parallel }{}^{final}(t)$ são as velocidades na direção da polarização, respectivamente no início e no final do interval Δt. No final deste intervalo, supomos que a direção de polarização muda, de $\overrightarrow{p}(t)$ para $\overrightarrow{p}(t+\mathrm{\Delta }t)$, e a velocidade inicial do intervalo de tempo subsequente é a projeção da velocidade $v_{\parallel }{}^{final}(t)\overrightarrow{p}(t)$ sobre o novo vetor polarização $\overrightarrow{p}(t+\mathrm{\Delta }t)$, isto é

(8) $$v_{\parallel }{}^{inicial}(t+\mathrm{\Delta }t)=v_{\parallel }{}^{final}(t)(\overrightarrow{p}(t)\cdot \overrightarrow{p}(t+\mathrm{\Delta }t)).$$

Na equação (8) supomos que a dinâmica da rede de actina é sujeita a um ruído que muda aleatoriamente a direção de polarização, como ditado pela equação (6). Supomos também que essa mudança de orientação reduz a velocidade paralela à polarização, já que $|\overrightarrow{p}(t)\cdot \overrightarrow{p}(t+\mathrm{\Delta }t)|\le 1$. Essa hipótese está de acordo com a proposição da correlação entre velocidade de migração e a organização do citoesqueleto, como proposto por Maiuri e colaboradores [²⁷].

Em conjunto, as equações (7) e (8) representam a dinâmica para a velocidade paralela à polarização, v_∥(t), como

(9) $$v_{\parallel }(t+\mathrm{\Delta }t)\overrightarrow{p}(t+\mathrm{\Delta }t)$$ $$=\left[(1-\gamma \mathrm{\Delta }t)v_{\parallel }(t)+{\displaystyle {\int }_t^{t+\mathrm{\Delta }t}}{\xi }_{\parallel }(t)\mathit{d}t\right]$$ $$\times (\overrightarrow{p}(t)\cdot \overrightarrow{p}(t+\mathrm{\Delta }t))\overrightarrow{p}(t+\mathrm{\Delta }t).$$

Um ponto crítico para o rigor da equação (9) diz respeito à dinâmica difusiva para θ(t), representada pela equação (4a, 4b). No entanto, a quantidade na equação (9) é $\overrightarrow{p}(t)\cdot \overrightarrow{p}(t+\mathrm{\Delta }t)=cos\mathrm{\Delta }\theta (t)\sim 1-\frac{{(\mathrm{\Delta }\theta )}^2}{2}$. Como ⟨(Δθ)²⟩∼Δt, já que a evolução de θ(t) segue um processo de Wiener, é válido supor na equação (9) que $\overrightarrow{p}(t)$ é constante quando Δt é pequeno. Mais detalhes sobre a validade desta equação estão na Ref. [³¹].

Supomos que a dinâmica na direção perpendicular à polarização é difusiva. Assim, a equação não pode levar em conta a velocidade, mas sim o deslocamento, que deve seguir um processo de Wiener:

(10) $$[r_{\perp }(t+\mathrm{\Delta }t)-r_{\perp }(t)]\overrightarrow{n}(t)=\overrightarrow{n}(t){\displaystyle {\int }_t^{t+\mathrm{\Delta }t}}{\xi }_{\perp }(t)\mathit{d}t,$$

onde $\overrightarrow{n}(t)=(\sin (\theta (t)),-\cos (\theta (t)))$ é um vetor unitário perpendicular a $\overrightarrow{p}(t)$.

ξ_∥(t), ξ_⊥(t), e β_⊥(t) são ruídos brancos Gaussianos (com diferentes unidades). ξ_∥(t) é independente dos outros dois, e nós tomamos como hipótese que ξ_⊥(t) e β_⊥(t) são relacionados: supomos que as flutuações da dinâmica da rede de actina no lamelipódio são responsáveis tanto pela mudança de orientação do eixo ré-frente como pelos deslocamentos estocásticos na direção perpendicular à polarização. Isto é, supomos que

(11) $${\xi }_{\perp }(t)=\sqrt{q}{\beta }_{\perp }(t),$$

com $\sqrt{q}$ dado em unidades de comprimento. As intensidades dos termos de ruído são dadas por meio de seus segundos momentos, como segue:

(12a) $$⟨{\xi }_{\parallel }(t)⟩=0,⟨{\xi }_{\parallel }(t){\xi }_{\parallel }(t\text{´})⟩=g\delta (t-t\text{´}),$$

(12b) $$⟨{\beta }_{\perp }(t)⟩=0,⟨{\beta }_{\perp }(t){\beta }_{\perp }(t\text{´})⟩=2k\delta (t-t\text{´}),$$

(12c) $$⟨{\xi }_{\perp }(t)⟩=0,⟨{\xi }_{\perp }(t){\xi }_{\perp }(t\text{´})⟩=2qk\delta (t-t\text{´}),$$

onde g, k, e qk têm unidades de [comprimento]²/tempo³, 1/tempo, e [comprimento]²/tempo, respectivamente.

O modelo pode ser resumido da seguinte maneira: consideramos uma partícula com dois graus de liberdade espaciais e um grau de liberdade interno, representando a polarização, que quebra a simetria espacial. A partícula segue uma dinâmica à la Langevin para a velocidade na direção da polarização, enquanto segue uma dinâmica de Wiener para o deslocamento na direção ortogonal. Há duas fontes independents de ruído: uma que age sobre a velocidade na direção de polarização e uma segunda fonte que resulta em um deslocamento difusivo na direção perpendicular à polarização.

Os resultados do modelo foram obtidos numericamente e analiticamente, usando integrais estocásticas. Detalhes das soluções são mostrados na Ref. [³¹]. O primeiro resultado diz respeito à relaxação da velocidade quadrática média da componente paralela à polarização, que tende a um valor estacionário:

(13) $$⟨v_{\parallel }{}^2(t)⟩={\displaystyle \frac{g}{2(\gamma +k)}}+\left(v_{\parallel }{}_0{}^2-{\displaystyle \frac{g}{2(\gamma +k)}}\right)\exp [-2(\gamma +k)t],$$

onde o limite assintótico representa o valor estacionário, definido como

(14) $$⟨v_{\parallel estacion\text{á}rio}^2⟩=\frac{g}{2(\gamma +k)}.$$

O tempo de relaxação R é obtido parte da equação (13):

(15) $$R={(\gamma +k)}^{-1}.$$

Os resultados numéricos podem ser comparados com a predição das equações (14) e (15) observando que o deslocamento total médio que ocorre em um determinado intervalo de tempo ε é a soma vetorial dos deslocamentos paralelo e ortogonal à polarização. Para o segundo caso, o valor médio do deslocamento quadrático é difusivo, portanto proporcional a ε, isto é,

(16) $$⟨{|r_{\perp }(t+\epsilon )-r_{\perp }(t)|}^2⟩=2kq\epsilon ,$$

como obtido em [³¹], de tal maneira que a podemos escrever $⟨\frac{{|\overrightarrow{r}(t+\epsilon )-\overrightarrow{r}(t)|}^2}{{\epsilon }^2}⟩=\frac{g}{2(\gamma +k)}+(v_{\parallel }^{2_0}-\frac{g}{2(\gamma +k)})\exp [-2(\gamma +k)t]+2kq/\epsilon$. A Figura 5 mostra esta quantidade como obtida a partir das trajetórias produzidas pela solução numérica, descontando $\frac{g}{2(\gamma +k)}$, mostrando a estabilização do valor 2kq/ε e a relaxação exponencial.

Figura 5
Gráficos semi-log de $⟨\frac{{|\overrightarrow{r}(t+\epsilon )-\overrightarrow{r}(t)|}^2}{{\epsilon }^2}⟩-\frac{g}{2(\gamma +k)}versust$, para q = 0.1, g = 10, γ = 1 e k como indicado. R obtidos como na equação (14). Símbolos correspondem à solução numérica para 10 sequências diferentes de números aleatórios. As linhas correspondem à solução apresentada pela equação (13).

A solução para o deslocamento quadrático médio pode ser escrita como

(17) $$\begin{array}{l}MSD=⟨|\overrightarrow{r}(t+\text{Δ}t)-\overrightarrow{r}(t){|}^2⟩\\ =\frac{g}{(\gamma +2k)(\gamma +k)}\left[\text{Δ}t-\frac{1}{\gamma +2k}(1-e^{-(\gamma +2k)\text{Δ}t})\right]+2qk\text{Δ}t,\end{array}$$

que pode ser reescrita na mesma forma que a Equação de Fürth modificada, equação (3), se identificarmos

(18a) $$D={\displaystyle \frac{g}{2(\gamma +2k)(\gamma +k)}},$$

(18b) $$P={\displaystyle \frac{1}{\gamma +2k}}$$

e

(18c) $$S=\frac{2qk(\gamma +2k)(\gamma +k)}{g+2qk(\gamma +2k)(\gamma +k)}.$$

Se, por um lado, a equação de Fürth modificada foi obtida como uma correção empírica, o modelo Anisotrópico tem como resultado analítico exato a mesma forma para o deslocamento quadrático médio. Isto significa que, por meio das equações (18a, 18b, 18c) podemos obter as quantidades macroscópicas P,D, e S, a partir dos parâmetros microscópicos do modelo.

Experimentos em migração celular consideram medidas de velocidade e de autocorrelação de velocidade. O modelo, no entanto, supõe que velocidade instantânea não é uma quantidade bem definida porque a componente na direção ortogonal à polarização diverge no limite de intervalo de tempo indo a zero. Por outro lado, as medidas experimentais de velocidade são sempre realizadas em intervalos de tempo finitos, tratando-se portanto de velocidades médias e não instantâneas. Em outras palavras, os resultados experimentais dependem do intervalo de tempo usado. No que segue avaliamos os efeitos do intervalo de tempo nas medidas de velocidade média e autocorrelação de velocidades.

Primeiro vamos usar a definição de velocidade instantânea em unidades naturais, $\overrightarrow{u}(\tau )$, como:

(19) $$\begin{array}{l}\overrightarrow{u}(\tau )=\underset{\delta \to 0}{\lim }\frac{\overrightarrow{\rho }(\tau +\delta )-\overrightarrow{\rho }(\tau )}{\delta }=\underset{\delta \to 0}{\lim }\frac{\text{Δ}{\overrightarrow{\rho }}_{\parallel }+\text{Δ}{\overrightarrow{\rho }}_{\perp }}{\delta }\\ =u_{\parallel }(t)\overrightarrow{p}(t)+\underset{\delta \to 0}{\lim }\frac{\text{Δ}{\rho }_{\perp }}{\delta }\overrightarrow{n}(t),\end{array}$$

onde $\mathrm{\Delta }{\overrightarrow{\rho }}_{\parallel }$ e $\mathrm{\Delta }{\overrightarrow{\rho }}_{\perp }$ são deslocamentos nas direções, respectivamente, paralela e ortogonal à polarização, medidos em unidades naturais, como definidas nas equações (4a, 4b). Quando k > 0 and q > 0, no limite em que δ→0, u_∥(τ) é bem comportada, mas a velocidade na direção orthogonal diverge, isto é , ${lim}_{\delta \to 0}\frac{\mathrm{\Delta }{\rho }_{\perp }}{\delta }$ vai para infinito, já que Δρ_⊥ segue um processo de Wiener e é proporcional a $\sqrt{\delta }$ . Por outro lado, em um experiment $\mathrm{\Delta }{\overrightarrow{\rho }}_{\parallel }$ e $\mathrm{\Delta }{\overrightarrow{\rho }}_{\perp }$ não são medidas separadamente. De fato, a re-análise dos dados de migração celular para 12 diferentes experimentos de 12 laboratórios apresentada por Thomas e colaboradores [¹⁹] mostra que as medidas do módulo para a velocidade estimada em intervalos de tempo finitos aumentam se o intervalo diminui.

Medidas de autocorrelação de velocidade (VACF, do inglês Velocity AutoCorrelation Function) podem trazer informação sobre a persistência dos movimentos. A função de autocorrelação de velocidades é definida de uma forma geral como

(20) $$VACF(\mathrm{\Delta }t)=⟨\overrightarrow{v}(t+\mathrm{\Delta }t)\cdot \overrightarrow{v}(t)⟩.$$

No entanto, se velocidade instantânea não está bem definida, é preciso redefinir como

(21) $$VACF(\mathrm{\Delta }t)=⟨v_{\parallel }(t)\overrightarrow{p}(t)\cdot v_{\parallel }(t+\mathrm{\Delta }t)\overrightarrow{p}(t+\mathrm{\Delta }t)⟩.$$

Neste caso, o resultado analítico fica [³¹]:

(22) $$VACF(\mathrm{\Delta }t)=⟨v_{\parallel estacion\text{á}rio}^2⟩e^{-\mathrm{\Delta }t/P}.$$

Por outro lado, nem sempre os resultados experimentais têm precisão suficiente para reproduzir o comportamento da equação (22). Há pelo menos dois efeitos possíveis, como discutimos no que segue.

Experimentos determinam a velocidade das células a partir de fotografias das suas localizações, o que implica em um intervalo de tempo finito. Seja δ este intervalo de tempo (em unidades naturais). Neste caso, pretende-se estimar a velocidade instantânea a partir da velocidade média $\overline{\overrightarrow{u}(\tau ,\delta )}$, definida como

(23) $$\overline{\overrightarrow{u}(\tau ,\delta )}={\displaystyle \frac{\overrightarrow{\rho }(\tau +\delta )-\overrightarrow{\rho }(\tau )}{\delta }}={\displaystyle \frac{\mathrm{\Delta }{\rho }_{\parallel }(\delta )}{\delta }}\overrightarrow{p}+{\displaystyle \frac{\mathrm{\Delta }{\rho }_{\perp }(\delta )}{\delta }}\overrightarrow{n}.$$

Podemos então definir uma função de autocorrelação para estas velocidades médias como

(24) $$\psi (\delta ,\mathrm{\Delta }\tau )=⟨\overline{\overrightarrow{u}(\tau ,\delta )}\cdot \overline{\overrightarrow{u}(\tau +\mathrm{\Delta }\tau ,\delta )}⟩,$$

que podem ser caculadas exatamente [³¹] e resultam

(25) $$\psi (\delta ,\mathrm{\Delta }\tau )={\displaystyle \frac{(\gamma +2k)}{\gamma }}{\displaystyle \frac{(1-e^{-\gamma \delta /(\gamma +2k)})(1-e^{-\delta })}{{\delta }^2}}\mathrm{}\times ⟨u_{\parallel sta}{}^2⟩e^{-\mathrm{\Delta }\tau }.$$

Quando δ > S, o deslocamento total da célula é medido no regime em que se observa um comportamento não difusivo, o que significa que Δρ_∥(δ) = u_∥δ≫Δρ_⊥. Neste caso, a velocidade média estimada é tal que $\overline{\overrightarrow{u}(\tau ,\delta )}\approx u_{\parallel }(\tau )\overrightarrow{p}(\tau )$ e estimar VACF usando $u_{\parallel }(\tau )\overrightarrow{p}(\tau )$ em vez da velocidade instantânea $\overrightarrow{u}(t)$ reproduz os resultados da equação (22).

Por outro lado, se o interval de tempo para a medida de velocidade é muito pequeno, de tal maneira que δ < S, o termo ortogonal à polarização domina e a estimativa para a velocidade torna-se $\overrightarrow{u}(\tau )\approx \frac{\mathrm{\Delta }{\rho }_{\perp }(\tau )}{\delta }\overrightarrow{n}(\tau )$, dando lugar a um valor VACF que vai a zero para valores decrescente de Δτ, já que Δρ_⊥ é governada por um processo de Wiener. Isto é de fato observado para experimentos, como mostrado por Thomas e colaboradores [¹⁹]. Estes desvios dependem de quão precisas são as medidas de velocidade, de tal maneira que a precisão da medida tem efeitos na curva experimental de ψ(δ,Δτ).

Considerando medidas de precisão infinita, resta ressaltar que mesmo para valores muito pequenos de δ, $\psi (\delta ,\mathrm{\Delta }\tau )\sim ⟨u_{\parallel sta}{}^2⟩e^{-\mathrm{\Delta }\tau }$, isto é, ψ(δ,Δτ) tende a VACF(Δτ). Para medidas de precisão finita, no entanto, ψ(δ,Δτ) diminui com Δτ se Δτ < S, devido à estimativa $⟨u_{\parallel sta}{}^2⟩$. A Figura 6 mostra o comportamento de ψ(δ,Δτ) com Δτ considerando 15, 2 ou 1 dígito signifcativo na medida das velocidades médias: pouca precisão implica estimativas da correlação que descrescem com Δτ.

Figura 6
Gráfico log-log em unidades naturais de ψ(δ,Δτ) versus Δτ para q = 0.1, g = 10, γ = 1, k = 0.04405(S = 0.001), usando δ = 0.001 e precisões diferentes para o cálculo da velocidade média. Quando a precisão é pequena, para a medida de velocidade, observamos que ψ(δ,Δτ) quando Δτ diminui.

Quando δ é infinitesimal, há ainda que garantir Δτ > δ, para que os intervalos de tempo [τ,τ + δ] e [τ + Δτ,τ + Δτ + δ] não se sobreponham. Uma vez que usamos estes intervalos para estimar $\overline{\overrightarrow{u}(\tau ,\delta )}$ e $\overline{\overrightarrow{u}(\tau +\mathrm{\Delta }\tau ,\delta )}$, quando Δτ < δ a sobreposição destes intervalos introduz correlações entre os deslocamentos. Isto acontece mesmo quando a precisão das medidas é alta. Estas interpretações a respeito das medidas de autocorrelação só foram possíveis de explicar à luz do modelo Anisotrópico. Interpretações equivocadas a respeito do comportamento observado da autocorrelação de velocidade são comuns na literatura, como realizadas nas referências [³², ³³].

Em resumo, a adição de uma nova variável interna, a polarização, possibilitou a definição de dois eixos com dinâmicas diferentes, quebrando a simetria espacial inerente ao problema proposto por Lagevin e resumido pelas equações (2). Na direção ortogonal à polarização a dinâmica é equivalente a um problema Browniano puro e, portanto, aparece como um termo de caráter difusivo nas equações para o MSD. Este comportamento descarta a possibilidade da definição de velocidade instantânea vetorial e colocando em cheque grandezas como a autocorrelação de velocidades. A solução apresentada por este modelo, separando as direções conforme o seu alinhamento com a polarização, propõe uma solução conceitual para este problema além de sugerir um protocolo experimental de medida da velocidade e de suas correlações.

Por outro lado, este modelo considera para uma das direções uma dinâmica Browniana pura, que foi o cerne da questão da discussão no início do século XX: como aplicar as leis de Newton para partículas para as quais não se pode medir velocidade instantânea. O ponto é que células não são partículas, mas corpos extensos, deformáveis e ativos, no sentido que obtêm energia cinética de uma maquinaria interna. Faz sentido que não haja um tempo finito entre duas deformações do lamelipódio, porque o lamelipódio sofre deformações simultâneas em locais diferentes.

Um modelo computacional que reproduza a cinética celular representa uma ferramenta matemática importante para a investigação de diferentes comportamentos. O modelo Anisotrópico, com a equação de Fürth modificada e as função de autocorrelação da velocidade dadas em unidades naturais fornece um protocolo para traduzir unidades das simulações (tipicamente pixels e passos de simulação) em unidades de laboratório e nos valores dos parâmetros micorscópicos, como veremos a seguir.

5. Modelo em CompuCell3D para Migração Celular

O modelo de simulação que discutimos nesta seção foi proposto por Fortuna e colaboradores [³⁴]. Está baseado no modelo de Potts Celular [³⁵, ³⁶, ³⁷] e foi realizado no ambiente Compucell3D [³⁸]. O modelo considera o espaço tridimensional dividido em uma rede cúbica. A cada sítio da rede são atribuídos dois índices. O primeiro índice indica se o sítio pertence a uma célula, ou se faz parte do meio ou do substrato. O segundo índice indica que tipo de compartimento celular aquele sítio compõe: núcleo, citoplasma ou lamelipódio. Nesta simulação, a célula é representada pelos domínios conexos de sítios com os índices correspondentes a núcleo, citoplasma e lamelipódio.

A dinâmica do modelo consiste em escolher um sítio e um dos seus vizinhos, e então tentar copiar os índices do sítio sobre os do vizinho, segundo um protocolo de Monte Carlo. Para obtar a probabilidade de que uma troca de índices seja aceita, calcula-se a variação da energia devido à cópia, considerando os seguintes termos:

(26) $$E=E_{interface}+E_{volume}+E_{F-actina},$$

onde E_interface é tal que cresce com as superfícies entre os vários componentes do sistema e, portanto, tenderá a minimizar as superfícies. Além disso, responde a uma hierarquia de valores tal que o espalhamento da célula sobre o substrato é favorecido (a célula ‘molha’ o substrato [³⁴]). E_volume é dada por

(27) $$E_{volume}={\lambda }_{cito}{\left(V_{cito}-V_{cito}^{alvo}\right)}^2+{\lambda }_{lamel}{\left(V_{lamel}-V_{lamel}^{alvo}\right)}^2+{\lambda }_{n\acute{u}cleo}{\left(V_{n\acute{u}cleo}-V_{n\acute{u}cleo}^{alvo}\right)}^2,$$

onde os subíndices ‘cito’, ‘lamel’ e ‘núcleo’ representam os componentes da célula (citoplasma, lamelipódio e núcleo) e têm o papel de garantir que cada componente da célula não possa variar demais seus volumes V_cito, V_lamel e V_núcleo, embora possibilite pequenas flutuações ao redor dos volumes alvo $V_{cito}^{alvo}$, $V_{lamel}^{alvo}$ e $V_{n\acute{u}cleo}^{alvo}$, que são tomados como parâmetros fixos das simulações. As constantes λ_cito, λ_lamel e ${\lambda }_{n\acute{u}cleo}$ regulam a intensidades destas flutuações. Finalmente, a variação de energia referente ao último termo de energia é calculada somente quando o sítio escolhido for um sítio de lamelipódio sobre o substrato e na borda com o meio, da seguinte forma

(28) $$\mathrm{\Delta }{E}_{F-actina}={\lambda }_{F-actina}[F(\overrightarrow{v})-F(\overrightarrow{r})],$$

onde $\overrightarrow{r}$ representa a localização do sítio escolhido (de lamelipódio) e $\overrightarrow{v}$, a posição do seu vizinho. $F(\overrightarrow{r})$ é um campo gerado em cada sítio de lamelipódio e que decai rapidamente quando o sítio não for lamelipódio, significando que $F(\overrightarrow{r})$ acaba por ter um perfil praticamente constante dentro do lamelipódio e que cai a zero quando fora. Assim, se na localização $\overrightarrow{v}$ está um sítio do meio, a diferença de energia é grande e negativa, favorecendo a cópia, mas se for um sítio de lamelipódio, ΔE_F−actina≈0. O subscrito F-actina neste termo refere-se aos filamentos de actina que formam a rede de actina polimerizada no lamelipódio.

A cópia dos índices do sítio escolhido será aceita se ΔE≤0 e, se ΔE > 0 ainda pode ser aceita, com probabilidade exp(−ΔE/T_B), onde T_B é um termo de flutuação à la Boltzmann. Mais detalhes sobre o código estão disponíveis na Ref. [³⁴].

O intuito das simulações é reproduzir a cinemática observada experimentalmente em migração de células únicas. Assim, as simulações iniciam com uma célula representada pelo núcleo rodeado de citoplasma, mas sem lamelipódio. Os sítios de citoplasma que tocam o substrato têm uma probabilidade finita de se transformar em lamelipódio enquanto o volume de lamelipódio for menor que seu volume alvo. A fração de volume alvo de lamelipódio, ϕ_l, em relação ao volume total da célula, é um dos parâmetros do modelo cuja variação foi investigada. Logo que o lamelipódio é formado, existe uma simetria radial no plano do substrato, o lamelipódio circunda a célula em todo o seu perímetro. A Figura 7 representa estas duas configurações, nos paineis A e B (visto de cima), A′ e B′ (vistos de lado e acima do plano do substrato) e A″ e B″ como um corte vertical, mostrando o núcleo em amarelo envolvido pelo citoplasma em cinza.

Figura 7
Início da migração para uma célula não polarizada, entrando em contacto com o substrato. Citoplasma em cinza, lamelipódio em verde e núcleo em amarelo. A: Célula não móvel na configuração incial, sem lamelipódio. B: célula não móvel, simétrica, com lamelipódio; C: competição entre lamelipódios incipientes, que acaba por levar à quebra de simetria; D: um dos lamelipódios acaba sendo selecionado e a célula torna-se capaz de movimento persistente. A célula transiciona entre C e D intermitentemente, mudando a direção do movimento em cada ciclo. Pontos de vista: A, B, C, D: vista superior; A′, B′, C′, D′: vista de um ponto acima e deslocada da célula; A″, B″, C″, D″: corte vertical. Figura da Ref. [³⁴].

No entanto, a dinâmica implícita pela equação (28) age como uma excitação local já que se há uma maior quantidade de lamelipódio, a probabilidade de se criar ali mais lamelipódio cresce. Por outro lado, o vínculo sobre o volume do lamelipódio representado pela equação (27), age como uma inibição global, diminuindo a criação de mais lamelipódio ali e em outros sítios. Este par – excitação local e inibição global – são ingredientes necessários para a quebra de simetria, como apontado por Turing [30]. Assim, a simetria quebra-se espontaneamente forma-se mais lamelipódio em um dos lados da célula, como mostrado nos paineis D, D′ e D″ da Figura 7. A formação de uma região com mais lamelipódio em um dos lados cria um eixo preferencial de cópias, de tal forma que a célula acaba por se mover nesta direção. Finalmente, locais diferentes podem apresentar lamelipódios incipientes (Figura 7, C, C′ e C″) e, eventualmente, pode acontecer de desestabilizar o lamelipódio dominante, mudando a orientação do movimento estocasticamente.

Para validarmos o modelo de simulação há que compará-lo com dados experimentais. Por meio do ajuste dos dados de simulação usando a equação de Fürth modificada, é possível repetir a análise dos dados experimentais mostradas anteriormente. Como as curvas após o ajuste podem ser colocadas na forma adimensional, dependente apenas do parâmetro S,a comparação de dados experimentais e resultados da simulação é direta. Fortuna e colaboradores [³⁴] simularam células com volume equivalente a esferas de raios $R_{c\acute{e}lula}=10$, 15 e 20 voxels, com frações volumétricas de lamelipódio variando de 0.05 a 0.30. Consideraram também diferentes intensidades para o campo de actina, representadas pelo valor de λ_F–actin. Os resultados para o MSD são mostrados na Figura 8 já ajustados e em unidades naturais. Os valores de S para estas simulações variam entre 2.59×10⁻⁵ até 6.25×10⁻³ [³⁴]. Nos insets das figuras, os dados são mostrados antes do ajuste.

Figura 8
MSD $(⟨{|\mathrm{\Delta }\overrightarrow{\rho }|}^2⟩)$ para a célula simulada com $R_{c\acute{e}lula}=15voxels$ em escala log-log e em unidades naturais. Médias sobre 5 rodadas de simulação. Os insets apresentam os mesmos dados em unidades de simulação: comprimento em múltiplos de raios celulares e tempo em passos de Monte Carlo. As curvas colapsam sobre o plot da Equação de Fürth modificada. As linhas tracejadas correspondem à equação de Fürth original. Adaptada da Ref. [³⁴].

As simulações possibilitam também a comparação dos resultados relacionados com a velocidade e autocorrelação de velocidade. A Figura 9 mostra que as autocorrelações de velocidade desviam-se do valor calculado pela equação (22), quando δ decresce. Este efeito é notável para δ < S. Estes resultados concordam com os resultados experimentais, mostrados na Figura 2, e com as previsões do modelo teórico, como mostrado na Figura 5.

Figura 9
O intervalo de tempo afeta δ usado para estimar a velocidade afeta o comportamento da função de autocorrelação de velocidade. Simulações $R_{c\acute{e}lula}=15$ voxels e ϕ_l = 0.10 e com dois valores para λ_F−actina = 150 e 200, que resultam valores de S fossem bem diferentes. As curvas ψ(δ,Δτ) como função de Δτ para diferentes valores δ mostram que se δ < S a autocorrelação cai se o intervalo de tempo é pequeno, reproduzindo os resultados experimentais e a previsão teórica.

Concluindo, este modelo de simulação

1. reproduz o mecanismo de excitação local e inibição global necessária para haver quebras espontâneas de simetria;

2. apresenta perda e reformação intermitente dos lamelipódios, típicos de migração de células únicas em substratos planos in vitro.

3. Reproduz a cinética observada experimentalmente nas curvas de deslocamento quadrático médio – MSD – nas funções de autocorrelação das velocidades, ψ_δ(Δτ).

Adicionalmente, este modelo computacional é suficientemente rápido para que possa simular situações mais complexas. Uma vez que o modelo foi validado pela comparação com dados experimentais, tem respaldo teórico e propõe protocolos de medida robustos, podemos utilizá-lo para investigar o efeito de campos químicos, substratos geometricamente complexos e migração coletiva de um conjunto de células.

6. Conclusões e Mensagens para Levar

Neste artigo apresentamos um problema que surgiu pela re-análise de dados experimentais de migração de célula única in vitro, sobre substratos planos. Os dados mostram um desvio do pardigma vigente para explicar este fenômeno, sob a forma de um regime difusivo em intervalos de tempo curtos, desafiando a hipótese que o modelo de Langevin, utilizado para descrever o movimento difusivo de partículas pontuais e inertes imersas em um fluido. Como vimos, os desvios mostrados nos dados experimentais não contradizem a aplicação do modelo de Langevin para partículas inertes e pontuais, mas a sua aplicação para corpos extensos, maleáveis e ativos.

Explicações anteriores buscaram explicar estes desvios por soluções simples demais (erros nas medidas de posição, o que representaria um termo constante adicionado à Equação de Fürth), ou complicados demais, propondo mecanismos adicionais que explicassem às vezes a perda de correlação ou ao aumento da correlação entre as velocidades. Neste último caso, aceitando a possibilidade destes mecanismos adicionais, os modelos poderiam ser considerados como complexos.

A re-análise de Thomas e coloaboradores deixou claro que havia problemas com o protocolo de medida de velocidade utilizado e, por consequência, com a definição da quantidade velocidade. A observação de que migração celular implica necessariamente em quebra de simetria espacial, pela formação de um eixo de polarização, juntamente com o fato do lamelipódio ser deformável e extenso, sugeriram a proposição do modelo Anisotrópico com uma direção onde o efeito de múltiplas deformações simultâneas e estocásticas em diferentes localizações do lamelipódio poderiam ser descritas pela deslocamento difusivo da célula na direção perpendicular à polarização e pela mudanças na direção de polarização.

O modelo de migração volta então a ser simples (embora requeira integração de processos estocásticos) se descrito pela Equação de Fürth modificada, já que se utiliza de três parâmetros macroscópicos (S, P e D). No entanto, este modelo se limita à cinética das células migrando isoladamente em substratos planos. As perguntas sobre efeitos de interação célula-célula, existência de obstáculos, diferentes topologias do substrato, campos químicos, ciclo celular estão todos fora do contexto do modelo. Se o nosso objetivo é investigar estes efeitos, precisamos de algum modelo mais complexo. O gás ideal termodinâmico é um modelo simples, mas se quisermos investigar os efeitos quânticos que sabemos existir, precisamos de mecânica quântica. O mesmo vale para fluidos: modelos para fluxos lineares são simples, mas se quisermos entender turbulência os modelos adequados podem ser muito complexos.

O estudo de migração celular volta-se agora para fenômenos que devem requerer modelos complexos onde surgem interações não lineares que limitam a resolução analítica dos modelos. Neste caso, estamos ainda tateando: não se sabe ainda quais os fatos que devem ser descartados e quais temos que conservar. Não sabemos quais medidas propor, no sentido se (1) serão robustas e (2) serão informativas. Por isso, não se tem certeza ainda de como projetar otimamente os experimentos. Por exemplo, um fenômeno importante no desevolvimento de seres multicelulares envolve migração coletiva. Para tanto diferentes experimentos têm sido propostos e diferentes medidas têm sido propostas. Mas variáveis de confusão estão sempre presentes e os mecanismos responsáveis nem sempre ficam claros.

O modelo em CompuCell3D, discutido acima e validado quantitativamente por experimentos e teoria, pode contribuir na determinação dos efeitos cinéticos de termos tais como interação de contacto, campos químicos, diferenças de adesão, etc. Simulações são mais rápidas, menos custosas e mais controláveis que experimentos e podem então sugerir medidas e comparações experimentais.

Finalmente, para enfatizar que a proposição de um modelo é uma tarefa para a qual não há receitas prontas (embora haja princípios a serem seguidos), termino parafraseando Tolstoi: os sistemas naturais bem modelados, são bem modelados da mesma maneira (fitando os dados experimentais). Os sistemas naturais (ainda) não bem modelados apresentam cada um a sua dificuldade.

Referências

• ^1.
  C.J. Weijer, Journal of Cell Science 122, 3215 (2009).
• ^2.
  P. Martin, Science 276, 75 (1997).
• ^3.
  E. Rognoni, A.O. Pisco, T. Hiratsuka, K.H Sipilä, J.M. Belmonte, S.A.Mobasseri, C. Philippeos, R. Dilão e F.M. Watt, Molecular Systems Biology14, e814 (2018).
• ^4.
  W.A. Muller, Trends in Immunology 24, 326 (2003).
• ^5.
  V. Mironov, V. Kasyanov e R.R. Markwald, Current Opinion in Biotechnology22, 667 (2011).
• ^6.
  C. Le Clainche e M.F. Carlier, Physiological Reviews 88, 489 (2008).
• ^7.
  F. Wang, Cold Spring Harbor Perspectives in Biology 1, a002980 (2009).
• ^8.
  L.P. Cramer, Nature Cell Biology 12, 628 (2010).
• ^9.
  A. Asnacius e O. Hamant, Trends in Cell biology 22, 584 (2012).
• ^10.
  A.T. Dawes e L. Edelstein-Keshet, Biophysical Journal 92, 744 (2007).
• ^11.
  N.W. Goehring e S.W. Grill, Trends in Cell Biology 72 (2013).
• ^12.
  R. Fürth, Zeitschrift für Physik 2, 244 (1920).
• ^13.
  P. Langevin, Comptes Rendus de l’Académie des Sciences de Paris 146, 530 (1908).
• ^14.
  A. Genthom, The European Physical Journal H 45, 49 (2020).
• ^15.
  A. Einstein, Annalen der Physik 17, 549 (1905).
• ^16.
  M. Smoluchowski, Bulletin International de l’Académie des Sciences de Cracovie1, 577 (1906).
• ^17.
  N. Wiener, Proceedings of the National Academy of Sciences 7, 294 (1921).
• ^18.
  F.A. Oliveira, R.M.S. Ferreira, L.C. Lapas e M.H. Vainstein, Front. Phys. 7 ,18 (2019).
• ^19.
  G.L. Thomas, I. Fortuna, G.C. Perrone, J.A. Glazier, J.M. Belmonte e R.M.C. Almeida,Physica A 550, 124493 (2020).
• ^20.
  P. Dieterich, R. Klages, R. Preuss e A. Schwab, PNAS 105, 459 (2008).
• ^21.
  H. Takagi, M.J. Sato, T. Yanagida e M. Ueda, PLoS One 3, e2648 (2008).
• ^22.
  P.H. Wu, A. Giri, S.X. Sun e D. Wirtz, PNAS 111, 3949 (2014).
• ^23.
  A.A. Potdar, J. Lu, J. Jeon, A.M. Weaver e P.T. Cummings, Annals of Biomedical Engineering 37, 230 (2009).
• ^24.
  L.E. Reichl, A modern course on Statistical Physics (WILEY-VCH, New Jersey, 2016),4^a ed.
• ^25.
  M. Boehm e E.G. Nabel, Circulation 103, 2879 (2001).
• ^26.
  J.L. Ross, PNAS 109, 5911 (2012).
• ^27.
  P. Maiuri, J.F. Rupprecht, S. Wieser, V. Ruprecht, O. Bénichou, N. Carpi, M.Coppey, S. De Beco, N. Gov, C.P. Heisenberg et al., Cell 161, 374 (2015).
• ^28.
  A.C. Callan-Jones e R. Voituriez, Current Opinion in Cell Biology 39, 12 (2016).
• ^29.
  M. Abercrombie, J.E.M. Heysman e S.M. Pegrum, Experimental Cell Research 62, 389 (1970).
• ^30.
  A.M. Turing, Phil. Trans. R. Soc. Lond. B 237, 37 (1952).
• ^31.
  R.M.C. Almeida, G.S. Giardini, M. Vainstein, J.A. Glazier e G.L. Thomas,arXiv:2002.00051 (2020).
• ^32.
  D. Selmeczi, S. Mosler, P.H. Hagedorn, N.B. Larsen e H. Flyvbjerg, Biophysical Journal 89, 912 (2005).
• ^33.
  L. Li, E.C. Cox e H. Flyvbjerg, Physical Biology 8, 046006 (2011).
• ^34.
  I. Fortuna, G.C. Perrone, M.S. Krug, E. Susin, J.M. Belmonte, G.L. Thomas, J.A.Glazier e R.M.C. Almeida, Biophysical Journal 118, 2801 (2020).
• ^35.
  F. Graner e J.A. Glazier, Physical Review Letters 69, 2013 (1992).
• ^36.
  J.A. Glazier e F. Graner, Physical Review E 47, 2128 (1993).
• ^37.
  P. Hogeweg, Biosystems 64, 97 (2002).
• ^38.
  M.H. Swat, G.L. Thomas, J.M. Belmonte, A. Shirinifard, D. Hmeljak e J.A. Glazier, Methods in Cell Biology 110, 32566 (2012).

Datas de Publicação

• Publicação nesta coleção
  05 Mar 2021
• Data do Fascículo
  2021

Histórico

• Recebido
  21 Set 2020
• Aceito
  17 Nov 2020