La caratterizzazione delle proprietà meccaniche della cellula è di grande importanza per stimarne lo stato biologico e la sua fisiologia. Tra le diverse tecniche di misura, tramite il microreometro è possibile di caratterizzarne le proprietà viscoelastiche.
In questo report si considera un esperimento di misura di un microreometro a biglia paramagnetica e il fitting di dati sperimentali per analizzare il comportamento di un modello meccanico di cellula a parametri concentrati.
Tale modello viene descritto dal parallelo di un corpo di Maxwell con una rigidezza elastica al quale viene aggiunto uno smorzatore in serie per ridurre la discrepanza con i dati sperimentali. Le costanti viscoelastiche stimate dal modello a parametri concentrati sono indicative delle proprietà dell’intera cellula e non dei singoli costituenti quali membrana, citoscheletro o citoplasma, per i quali sarebbe richiesto un modello differente.
Il seguente articolo è un estratto di un progetto realizzato per il corso di Modellazione e Simulazione di Sistemi Fisiologici. Tenuto per Ingegneria Medica presso l’Università degli Studi di Roma Tor Vergata.
Introduzione
Per caratterizzare meccanicamente la cellula e le sue prprietà viscoelastiche è necessario avere un modello descrittivo del suo comportamento. Viene quindi considerato un modello a parametri concentrati costituito da uno smorzatore in serie ad un ramo dato dal paralello di un corpo di Maxwell e una rigidezza elastica. Si considerano quindi dei parametri materiali noti dalla letteratura tramite i quali viene analizzata la soluzione, analitica e numerica, a differenti tipologie di ingresso.
Background
Le proprietà e la fisiologia cellulare sono il principale motivo di indagine e ricerca. Tra queste le proprietà meccaniche sono utili biomarker [1].
La deformabilità e le proprietà viscoeleastiche sono utili indicatori di cambiamenti citoscheletrici e dello stato cellulare. Ne sono esempi l’aumento di rigidezza nei globuli rossi affetti da sferocitosi [2, 3] e l’aumento di deformabilità nella cellule cancerogene [4]. Tali proprietà meccaniche possono essere misurate con differenti tecniche quali microscopio a forza atomica, trappole ottiche, microreometro a biglia magnetica e l’aspirazione con micropipette [5].
In questo report viene posta l’attenzione su un microreometro a biglie magnetiche e il modello di cellula per stimarne le proprietà viscoelastiche.
Microreometro a biglia magnetica
Il microreometro è un dispositivo utilizzato nella micro reologia per caratterizzare le proprietà reologiche di un mezzo quali la viscosità o le tracce di flusso. Tali dispositivi si misura si dividono in attivi e passivi a seconda che sfruttano l’energia termica o forze esterne applicate per movimentare la particella.
In questo report viene analizzato un microreometro a biglie magnetiche. Tali biglie paramagnetiche di 4.5 𝜇𝑚 di diametro sono legate al citoscheletro cellulare tramite la connessione tra le integrine di membrana e la fibronectina di cui sono ricoperte le biglie. L’integrina di membrana è connessa direttamente al citosceheltro quindi la misura dello spostamento della biglia fornisce informazioni sulla meccanica cellulare.
Magnetizzazione
La presenza di un campo magnetico controllato applica una forza sulla biglia e misurandone lo spostamento è possibile caratterizzarne alcune proprietà viscoelastiche [6] fissando un modello viscoelastico di cellula. Il comportamento osservato è rappresentato in fig. 1b ed è possibile distinguere tre fasi principali. Come viene applicata la forza è c’è un primo spostamento della biglia (fase I) identificato come risposta elastica veloce. La prima fase è seguita dal regime di rilassamento (fase II) e infine la regime di flusso, dove si verifica un progressivo scivolamento.
Considerando la seguente descrizione del modello è possibile stimarne sperimentalmente tre parametri quali l’ampiezza della prima fase 1/𝑘1+𝑘0 , la costante di tempo dell’andamento esponenziale della fase II pari a 𝜏 = 𝛾1(𝑘0+𝑘1)/𝑘0 𝑘1 e il parametro 𝛾0 identificativo dell’ampiezza della fase III. Tali valori, stimati da Bausch et al. [6], verranno utilizzati per implementare una soluzione numerica del modello all’interno dell’ambiente Matlab.
Modello meccanico
Tali curve sperimentali sono analizzate sulla base di un modello meccanico descritto da una serie di uno smorzatore (a seguire indicato come Dashpot) e di un corpo Zener, noto anche come corpo Kelvin, descritto da un parallelo di una ramo con una rigidezza elastica e un secondo ramo formato dalla serie di uno smorzato e una rigidezza. Tale modello è rappresentato in fig. 1a. La descrizione meccanica segue dal fatto che lo spostamento risultante sarà proprio la somma del dashpot 𝛾0 (D) e del corpo Zener (Z).
\begin{equation}
x(t)=x_{D}(t)+x_{z}(t)
\end{equation}E quindi derivando vale ancora:
\begin{equation}
\dot{x}(t)=\dot{x}_{D}(t)+\dot{x}_{z}(t)
\end{equation}Chiaramente la forza applicata sarà la stessa su entrambi.
Per (D) è nota la relazione dello smorzatore da cui possiamo ricavare la velocità come:
\begin{equation}
\dot{x}_{D}=\frac{F}{\gamma_{0}}
\end{equation}Per il corpo Zener la forza si distribuisce sul ramo elastico e sul ramo dove è presente un corpo Maxwell (molla e smorzatore).
\begin{equation}
F(t)=F_{\text {Maxwell }}(t)+F_{\text {Ko }}(t)
\end{equation}Dove la forza applicata sul corpo di Maxwell è proprio la differenza tra quella totale applicata e quella del ramo elastico 𝑘0𝑥𝑍(𝑡). Ovvero per il corpo Maxwell, derivando lo spostamento, si ottiene:
\begin{equation}
\dot{x}_{Z}(t)=\frac{\dot{F}}{k_{1}}+\frac{F}{\eta_{0}}
\end{equation}E quindi complessivamente per il corpo Zener vale:
\begin{equation}
\dot{x}_{Z}(t)=\frac{\dot{F}(t)-k_{0} \dot{x}(t)}{k_{1}}+\frac{F(t)-k_{0} x(t)}{\eta_{0}}
\end{equation}In conclusione, lo spostamento sarà dato dalla somma dei due spostamenti che possono essere ricavati risolvendo le due equazioni differenziali per un dato ingresso e condizione iniziale.
Risultati
Tale modello viene quindi analizzato all’interno dell’ambiente numerico Matlab confrontando la soluzione numerica e analitica per differenti tipologie di ingressi.
Ingresso a gradino
Per un ingresso a gradino, di ampiezza ̄ 𝐹 otterremo per il dashpot l’equazione differenziale:
\begin{equation}
\left\{\begin{array}{l}
\dot{x}_{D}=\frac{\bar{F}}{\gamma_{o}} \quad t>0 \\
x_{D}(0)=0
\end{array}\right.
\end{equation}Di soluzione analitica lineare nel tempo:
\begin{equation}
x_{D}(t)=\frac{\bar{F}}{\gamma_{0}} t
\end{equation}Tale risposta lineare è rappresentante della fase lineare osservata nelle misure sperimentali. Per il corpo Zener vale:
\begin{equation}
\left\{\begin{array}{l}
\dot{x}_{z}(t)\left(1+\frac{k_{o}}{k_{1}}\right)=-\frac{k_{o}}{\gamma_{1}} x_{z}+\frac{\bar{F}}{\gamma_{1}} \quad t>0 \\
x_{z(0)}=\frac{\bar{F}}{k_{1}+k_{o}}
\end{array}\right.
\end{equation}Dove la condizione iniziale deriva dal fatto da considerare lo spostamento sullo smorzatore nullo all’instante iniziale e allora lo spostamento complessivo è legato alle sole rigidezze elastiche.
function dydt=zener_displacement(t,y,parameters,tspan,flag)
k_0=parameters(1);
k_1=parameters(2);
gamma_0=parameters(3);
gamma_1=parameters(4);
F_bar=parameters(5);
switch flag
case ’square’
[F, dF]=force(t,tspan);
dydt=((F-k_0*y)/gamma_1 + (dF/k_1))*(1+k_0/k_1);
case ’step’
dydt=(-(k_0/gamma_1)*y+(F_bar/gamma_1))/(1+(k_0/k_1));
case ’harmonic’
omega=parameters(7);
F=F_bar*(1+sin(omega*t));
dF=F_bar*omega*cos(omega*t);
dydt=((F-k_0*y)/gamma_1 + (dF/k_1))*(1+k_0/k_1);
end
endfunction dydt=dashpot_displacement(t,y,parameters,tspan,flag) gamma_0=parameters(3);
gamma_1=parameters(4);
F_bar=parameters(5);
switch flag
case ’square’
[F, dF]=force(t,tspan);
dydt=F/gamma_0;
case ’step’
dydt=F_bar/gamma_0;
case ’harmonic’
omega=parameters(7);
F=F_bar*(1+sin(omega.*t));
dF=F_bar*omega*cos(omega*t);
dydt=F/gamma_0;
end
endConsiderando il tempo di rilassamento:
\begin{equation}
\tau=\gamma_{1} \frac{k_{1}+k_{1}}{k_{0} k_{1}}
\end{equation}Si può riscrivere il problema come:
\begin{equation}
\left\{\begin{array}{l}
\dot{x}_{z}(t)+\frac{1}{\tau} x_{Z}(t)=\frac{\bar{F}}{k_{0} \tau} \quad t>0 \\
x_{z(0)}=\frac{\bar{F}}{k_{1}+k_{o}}
\end{array}\right.
\end{equation}Ovvero, tale problema di Cauchy presenta una soluzione del tipo:
\begin{equation}
x_{z}(t)=e-{ }^{A(t)}\left[x_{0}+\int_{t_{0}}^{t} g(x) e^{A(t)}\right]
\end{equation}Dove 𝐴(𝑡) è il termine legato alla soluzione dell’equazione omogenea, ovvero l’integrale dell’esponenziale con esponente il coefficiente che moltiplica 𝑥𝑍(𝑡), mentre il termine 𝑔(𝑥) è legato all’equazione non omogenea. Allora si ottiene:
\begin{equation}
x_{z}(t)=e-\frac{t}{\tau}\left[\frac{F}{k_{0}+k_{1}}+\frac{F}{k_{0}} e^{\frac{t}{\tau}}-\frac{F}{k_{0}}\right]
\end{equation}Che può essere riscritta come:
\begin{equation}
x_{Z}(t)=\frac{\bar{F}}{k_{0}}\left(1-\frac{k_{1}}{k_{0}+k_{1}} e^{-\frac{t}{\tau}}\right)
\end{equation}Ovvero un’esponenziale che progressivamente tenderà a ̄ 𝐹/𝑘0.
Inoltre, lo spostamento è raffigurato in fig. 3a e fig. 9. Tale soluzione viene graficata in fig. 3a e confrontata con la soluzione numerica ottenuta in Matlab.
Soluzione numerica
La soluzione numerica viene ottenuta in Matlab tramite 𝚘𝚍𝚎𝟷𝟻𝚜, un solutore proprietario per equazioni differenziali stiff.
La routine principale (main.m) permette di selezionare tramite la variabile 𝚏𝚕𝚊𝚐 il tipo di ingresso per la forza, scegliendo tra l’ingresso a gradino, un’onda quadra e una sinusoide, descritte più avanti. Successivamente la routine setta un intervallo temporale appropriato e richiama 𝚘𝚍𝚎𝟷𝟻𝚜() passandogli le funzioni dove sono definiti gli spostamenti (fig. 5) e poi viene graficata la loro somma.
Ingresso ad onda quadra
Viene poi testato un ingresso ad onda quadra come utilizzato in [6], con un periodo di 5s, un duty cicle del 50% e un’ampiezza di 2000 pN, testata su un intervallo di 20s. Tale ingresso viene implementato tramite la routine force.m, descritta in fig. 8 avendo l’accortezza di ammorbidire il gradino per poterne fare la derivata numerica anche delle discontinuità e facendo in modo che 𝚘𝚍𝚎𝟷𝟻𝚜 possa richiederne il valore a qualsiasi instante, e non solo sui campione dove viene definita, tramite un’interpolazione.
Tale segnale è presente in fig. 7. I risultati sono presenti in fig. 3b e dalla scomposizione dei due segnali (fig. 9) è possibile osservare come il dashpot presenta una componente lineare a tratti alla quale viene aggiunto il comportamento oscillatore del corpo Zener.
Ingresso sinusoidale
L’effetto oscillatorio è chiaramente molto presente se sollecitiamo il modello con un ingresso sinusoidale del tipo:
\begin{equation}
F(t)=\bar{F}(1+\sin (\omega t))
\end{equation}La soluzione analitica presenta sia una componente armonica che una esponenziale:
\begin{equation}
\begin{aligned}
x(t)=\frac{F_{0}}{k_{0}} & {\left[1-\frac{k_{1}}{k_{0}+k_{1}} e^{-t / \tau}+\right.} \\
&+\frac{\tau \omega\left(1-\frac{\gamma_{1}}{k_{1} \tau}\right)\left(e^{-t / \tau}-\cos \omega t\right)}{1+(\tau \omega)^{2}}+\\
&\left.+\frac{\left.\left(1+\frac{\gamma_{1} \tau \omega^{2}}{k_{1}}\right) \sin \omega t\right]}{1+(\tau \omega)^{2}}\right]
\end{aligned}
\end{equation}Tale soluzione viene confrontata con la soluzione numerica in fig. 6. Dai due spostamenti separati (fig. 9) è possibile vede come entrambi hanno una componente armonica e mentre lo spostamento del corpo Zener è limitato il dashpot presenta una crescita lineare inviluppato da una sinusoide.
function [F, dF]=force(t_curr,t)
T_cyc =5;
N_cyc =4;
F=2000;
freq=1/T_cyc;
squareWave=F/2*(square(2*pi*freq*t,50)+1);
squareWave(end)=0;
squareWave_smooth=smoothdata(squareWave ,’gaussian’,(size(t)/N_cyc)*0.05);
derivate=[0 diff(squareWave_smooth)];
F=interp1(t,squareWave_smooth ,t_curr);
dF=interp1(t,derivate ,t_curr);
end
Conclusioni
Lo studio e la classificazione del comportamento meccanico di tali modelli permette di analizzarne i coefficienti e poter fittare dei dati sperimentali per ricavare le proprietà meccaniche cellulari. Avvalersi di un solutore numerico permette di analizzare il modello meccanico anche per ingressi complessi e fortemente variabili nel tempo.
Le costanti viscoelastiche stimate dal modello a parametri concentrati sono indicative delle proprietà dell’intera cellula, ma tali proprietà sono determinate in parte dalle proprietà della membrana cellulare e in parte dalle proprietà del citoplasma e dalla struttura del citoscheletro. Per avere maggiori informazioni sulle singole componenti cellulari è richiesto un modello più dettagliato.
Disponibilità del codice
Il codice è disponibile alla repository GitHub del progetto.
Riferimenti
- [1] Dino Di Carlo. “A Mechanical Biomarker of Cell State in Medicine”. en. In: SLAS Technology 17.1 (Feb. 2012), pp. 32–42. ISSN: 24726303. DOI: 10 . 1177 / 2211068211431630.
- [2] Gabriel Y.H. Lee and Chwee T. Lim. “Biomechanics approaches to studying human diseases”. en. In: Trends in Biotechnology 25.3 (Mar. 2007), pp. 111–118. ISSN: 01677799.
- [3] S. Suresh et al. “Connections between single-cell biomechanics and human disease states: gastrointestinal cancer and malaria”. en. In: Acta Biomaterialia 1.1 (Jan. 2005), pp. 15–30. ISSN: 17427061.
- [4] Claudia Tanja Mierke. “Viscoelasticity Acts as a Marker for Tumor Extracellular Matrix Characteristics”. In: Frontiers in Cell and Developmental Biology 9 (Dec. 2021), p. 785138.
- [5] C Ross Ethier and Craig A Simmons. Introductory Biomechanics: From Cells to Organisms. en.
- [6] Andreas R. Bausch et al. “Local Measurements of Viscoelastic Parameters of Adherent Cell Surfaces by Magnetic Bead Microrheometry”. en. In: Biophysical Journal 75.4 (Oct. 1998), pp. 2038–2049. ISSN: 00063495.
