Diffusione e FEA per sfere a rilascio di farmaco

La modellazione matematica dei meccanismi di diffusione per il rilascio di farmaci può essere molto utile per accelerare lo sviluppo del prodotto e per comprendere meglio i meccanismi che controllano il rilascio di farmaci da sistemi di somministrazione avanzati.

L’integrazione dell’ambiente numerico di Comsol con la programmazione in Matlab permette di sviluppare modelli risolutivi personalizzati e facilmente integrabili. Tali modelli permettono agevolmente l’analisi di diffusione per sistemi a rilascio di farmaci con geometrie e proprietà costitutive desiderate.

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

I sistemi a rilascio di farmaco sono sempre più diffusi e la simulazione in silico può essere notevolmente di aiuto nella scelta del farmaco e delle procedure di produzione al fine di fornire un profilo di rilascio specifico.

Grazie ai significativi progressi nella tecnologia dell’informazione, la modellazione matematica della somministrazione di farmaci è un campo di importanza accademica e industriale in costante aumento con un enorme potenziale futuro.

Ci si può aspettare che l’ottimizzazione in silico di nuovi sistemi di somministrazione di farmaci aumenti in modo significativo l’accuratezza e la facilità di applicazione. Inoltre, è probabile che le simulazioni al computer diventino parte integrante della futura ricerca e sviluppo nella tecnologia farmaceutica.

Risulta quindi fondamentale avere la possibilità di effetuare analisi per la composizione, geometria, dimensioni e la procedura di preparazione dei vari tipi di sistemi di somministrazione, tenendo conto della via di somministrazione desiderata, della dose del farmaco e del profilo di rilascio.

Background

Il rilascio controllato di sostanze trova origine intorno al 1950 con lo sviluppo di sistemi polimerici capaci di effettuare un rilascio programmato [1].

Sistemi a rilascio di farmaco

Tale controllo può essere applicato per diversi ambiti quali cosmetica, industria farmaceutica, agricoltura e scienze alimentari.

In particolare, concentrandosi su sistemi per il rilascio di farmaci questi permettono di avere un rilascio graduale e continuo nel tempo garantendone la dose all’interno dei range ottimali. Tali sistemi possono essere sviluppati in modo tale da ottimizzare la cinetica del rilascio garantendo la massima efficacia della terapia farmacologica.

Esistono diverse tipologie di sistemi per il rilascio controllato basati su differenti fenomeni quali diffusione, dissoluzione, scambio ionico e controllo con mediatori chimici [2].

Dunque, in questo report viene posta l’attenzione su sistemi basati sulla diffusione. Inoltre, tali sistemi si basano sull’assunzione che il trasporto di massa limita la velocità di rilascio del farmaco, le condizioni modellate sono tali per tutta la durata del rilascio, il dispositivo non si sgonfia eccessivamente durante il rilascio e il coefficiente di diffusione rimane costante nello spazio e nel tempo [3].

Diffusione-convezione

La diffusione è la migrazione di una sostanza da una regione ad alta concentrazione ad una a concentrazione minore. I farmaci sono intrappolati all’interno di capsule inerti insolubili in acqua costituite da membrane polimeriche (approccio a sistema di riserva) oppure all’interno di una matrice polimerica (approccio monolitico). Pertanto in entrambi i casi il farmaco diffonde dalla sfera carica alla soluzione.

Pertanto i meccanismi di diffusione sono regolati, nello spazio e nel tempo, dalla legge di Fick. Tale legge lega la densità di corrente di soluto al gradiente di concentrazione, per il tramite del coefficiente di diffusione:

\begin{equation}
\mathbf{j}_{d}=-D \nabla c
\end{equation}

A cui si aggiunge un contributo convettivo:

\begin{equation}
\mathbf{j}_{c}=c \mathbf{q}
\end{equation}

Quindi focalizzandosi su un volume di controllo si può considerare il bilancio di massa. La massa sarà l’integrale spaziale della concentrazione:

\begin{equation}
m_{p}(t)=\int_{P} c(\mathbf{x}, t) d v
\end{equation}

Ovvero tale approccio consente di portare la derivata all’interno dell’integrale, tramite il teorema di Reynolds.

\begin{equation}
\frac{d m(t)}{d t}=\int_{p} \frac{\partial c(\mathbf{x}, t)}{\partial t} d v=-\int_{\partial P} \mathbf{j} \cdot \hat{n} d a+\int_{P} \mathbf{r} d v
\end{equation}

Pertanto la variazione della massa sarà legata ai contributi entranti dal bordo più un eventuale contributo di sorgente. Ovvero, applicando il teorema della divergenza, la conservazione della massa può essere riscritta come:

\begin{equation}
\int_{P}\left(\frac{\partial c}{d t}+\operatorname{div} \mathbf{j}\right) d v=\int_{P} \mathbf{r} d v
\end{equation}

Localizzazione

Quindi, sfruttando l’arbitrarietà dal volume di controllo, questo può essere localizzato:

\begin{equation}
\frac{\partial c}{\partial t}+\operatorname{div} \mathbf{j}=\mathbf{r}
\end{equation}

E, complessivamente, per un fluido incomprimibile, si ottiene un problema differenziale:

\begin{equation}
\left\{\begin{array}{l}
\frac{\partial c}{\partial t}-D \Delta_{2} c+q \cdot \nabla \mathbf{c}=\mathbf{r} \\
\text { +B. C. } \\
\text { +I. C. }
\end{array}\right.
\end{equation}

Devono dunque essere considerate anche le condizioni iniziali, ovvero la distribuzione delle specie nel sistema, e le condizioni al bordo. Allora tale equazione alle derivate parziali può essere risolta agevolmente con un solutore agli elementi finiti.

Inoltre, per la descrizione di un eventuale rivestimento del volume si potrebbe considerare un guscio tridimensionale sufficientemente sottile con un suo coefficiente di diffusione. Tuttavia, è più conveniente considerare un’interfaccia, ovvero una superficie 2D.

Per la descrizione vale un’ammettenza di coating per la quale vale:

\begin{equation}
\mathbf{j} \cdot(-\hat{n})=h[[c]]
\end{equation}

Spesso, tale valore si considera essere:

\begin{equation}
h=\frac{D_{\text {barriera }}}{\text { spessore }}
\end{equation}

Alla quale si aggiunge l’equazione di continuità:

\begin{equation}
[[\mathbf{j} \cdot(-\hat{n})]]=0
\end{equation}

Comsol Multiphysics

Tale modello di diffusione-convezione viene implementato all’interno dell’ambiente COMSOL Multiphysics [4].

Tale ambiente permette di effettuare simulazioni avvalendosi di metodi numerici avanzati e sfruttando l’interfaccia grafica tramite la quale interagire unendo fenomeni fisici differenti all’interno della stessa analisi.

Viene largamente utilizzato in campo scientifico e ingegneristico per simulare e progettare dispositivi e processi in molti scenari applicativi. Oltre all’applicazione di creazione modelli viene sfruttata l’integrazione con l’ambiente Matlab trasferendo l’impostazione e la soluzione del modello sulla programmazione Matlab.

Risultati

Per avere un modello agevole da impostare ed eventualmente anche da modificare, viene impostata una struttura a blocchi, ciascuno definito da una subroutine. Lo schema a blocchi è presente in fig. 3.

Il codice parte caricando i pacchetti Comsol e attivando il percorso verso la cartella contenente le subroutines. Dopo la prima fase di inizializzazione vengono impostati i parametri della simulazione, presenti in Tab. 1, caricati tramitemydata.m. Tali parametri vengono anche ricaricati in Comsol tramite set_parameters(). Questo, come altre accortezze nel seguito, permettono di riaprire tale file in Comsol e mantenere tutti i valori precedentemente salvati all’interno di Matlab.

Una volta caricati i parametri si può passare ad impostare la geometria. La routine create_geometry() si occupa di creare il modello geometrico all’interno di Comsol. Vengono create le N geometrie sferiche (𝑁 = 4 in questa analisi) e la geometria cilindrica considerando i parametri (posizione e dimensioni) precedentemente citati.

Tale routine permette anche la selezione dei domini e dei bordi delle geometrie. Ciascun elemento geometrico creato definisce un dominio interno e alcuni elementi di bordo tramite un ID numerico che risulta necessario per applicare condizioni al bordo o condizioni iniziali.

Modello numerico

In particolare, la tecnica per identificarli sfrutta la costruzione di un bounding box intorno i singoli elementi e la richiesta, al modello Comsol, dei domini (o bordi) al suo interno.

Infine, tale routine restituisce il modello aggiornato con le geometrie e gli ID dei domini e dei bordi di sfere e cilindro.

Successivamente viene aggiunta la fisica tramite la routine create_physics(). Vengono aggiunte le proprietà di diffusione al cilindro e alle sfere, tramite il coefficiente di diffusione, considerato isotropo, e le condizioni iniziali sulle sfere.

Dunque viene meshata la geometria con più o meno accuratezza, scelta tramite un ID numerico da 1 (Extremerly Fine) a 0 (Extremely Coarse). Viene anche implementata una routine ad hoc per convertire la flag testuale nell’ID numerico.

Dipendenza temporale

Si crea poi lo studio numerico, si impostano il tempo di inizio e fine e la discretizzazione temporale per poi risolvere il modello. La routine risolutiva crea uno studio tempo dipendente e richiama i diversi solutori numerici proprietari di Comsol.

Infine, il modello viene passato alla routine per il post pocessing. La concentrazione viene integrata nel volume per ottenere la massa nel tempo, il cui andamento è riportato in fig. 6a. Viene anche letta la concentrazione in direzione radiale nelle sfere. Per farlo vengono create delle cutline ad hoc che tagliano le sfere dal centro al bordo. Successivamente vengono estratti i dati esattamente sulle cutline e sono riportati in fig. 6b.

Analisi di convergenza

Viene poi effettuata un’analisi di sensibilità al magliaggio andando a ripetere, all’interno di un loop for, l’analisi variando il grado di rifinimento della mesh. Vengo quindi raccolti i dati per ogni simulazione quali il numero di gradi di libertà e viene considerato come valore di interesse la massa finale residua nella sfere all’ultimo instante temporale. L’errore viene quindi calcolato rispetto il valore assunto veritiero, ovvero il dato della simulazione con mesh il più possibile fitta.

Dai cui risultati, in fig. 5a e in Tab. 2, si sceglie la mesh con ID numerico pari a 2. Si sceglie quindi una discretizzazione Extra Fine con 154456 gradi di libertà (fig. 5b).

Ovvero tale soluzione presenta un errore inferiore al 5% rispetto la soluzione considerata vera, ovvero la soluzione con la discretizzazione massima possibile, ma con un costo computazionale 5 volte inferiore.

Evoluzione della massa

Dunque, una volta risolta la simulazione numerica è possibile analizzare la variazione di diverse quantità.

A seguire, viene analizzato l’andamento della massa all’interno del soluto che progressivamente aumenta in funzione del rilascio di farmaco da parte delle sfere. L’andamento è riportato in fig. 6a. Allora la stessa quantità viene poi letta anche all’interno delle quattro sfere dove progressivamente diminuisce in modo simile tra loro (fig. 6b).

Dunque vengono letti anche i valori della concentrazione, in direzione radiale, all’interno delle quattro sfere. L’andamento risulta simile tra loro e viene riportato solo l’andamento per la sfera numero 4, in fig. 4. Vengono poi riportate anche le mappe colori di superficie con la concentrazione all’instante finale e finale, in fig. 7.

Coating superficiale

Viene poi aggiunta la possibile di inserire un coating superficiale, ovvero un’interfaccia sulla sfera che ne altera le proprietà di diffusione.

Tale peculiarità, i cui parametri sono presenti in Tab. 1, viene passato direttamente routine per la creazione della fisica indicando il numero delle sfere (da 1 a 4, rispetto la fig. 8a) rispetto il quale sono state inizializzate.

Successivamente la routine inserisce le proprietà di coating utilizzando l’opzione di barriera sottile e fornendo a Comsol i relativi parametri. All’utente viene poi restituito un feedback visivo tale da indicare le sfere sulle quale viene applicato il coating (fig. 8b).

Effetti del rivestimento

Tale rivestimento rallenta notevolmente la diffusione e ciò risulta evidente sia dalla massa totale nel cilindro (fig. 9a) che dalla massa residua nelle sfere (fig. 9b). Ovvero la massa nelle sfere progressivamente si riduce ma lo fa con velocità diverse.

Inoltre, analizzando l’andamento della concentrazione nelle sfere, in fig. 10 si vede come la concentrazione nelle sfere con il coating rimane più alta al passare del tempo.

Estensione temporale

Analizzando la fig. 9b è evidente come il problema risulta temporalmente incompleto per cui viene estesa l’analisi da 4 a 40 secondi mostrando come la massa nelle due sfere tende ad equilibrarsi allo stesso valore, seppur con velocità diverse. Analogamente la massa nel cilindro tende ad una situazione di equilibrio (fig. 11).

In appendice sono riportate le routine per l’estrazione numerica dei risultati dal modello Comsol.

Conclusioni

L’integrazione dell’ambiente numerico di Comsol con la programmazione in Matlab permette di sviluppare modelli risolutivi personalizzati e facilmente integrabili. Tali modelli permettono agevolmente l’analisi di diffusione per sistemi a rilascio di farmaci con geometrie e proprietà costitutive desiderate.

Lo studio del rilascio di farmaco e della sua diffusione permette di analizzarne la variazione nel tempo e le sue concentrazioni nella soluzione principale.

Riferimenti

• [1] Cong Truc Huynh and Doo Sung Lee. “Controlled Re- lease”. en. In: Encyclopedia of Polymeric Nanoma- terials. Ed. by Shiro Kobayashi and Klaus Müllen. Berlin, Heidelberg: Springer Berlin Heidelberg, 2015, pp. 439–449. ISBN: 978-3-642-29647-5 978-3-642- 29648-2. DOI: 10.1007/978-3-642-29648-2_314. URL: http://link.springer.com/10.1007/978- 3-642-29648-2_314 
• [2] Ronald A. Siegel and Michael J. Rathbone. “Overview of Controlled Release Mechanisms”. en. In: Fundamen- tals and Applications of Controlled Release Drug De- livery. Ed. by Juergen Siepmann, Ronald A. Siegel, and Michael J. Rathbone. Boston, MA: Springer US, 2012, pp. 19–43. ISBN: 978-1-4614-0880-2 978-1- 4614-0881-9. DOI: 10.1007/978-1-4614-0881- 9_2. URL: http://link.springer.com/10.1007/ 978-1-4614-0881-9_2 
• [3] Juergen Siepmann and Florence Siepmann. “Modeling of diffusion controlled drug delivery”. en. In: Journal of Controlled Release 161.2 (July 2012), pp. 351–362. ISSN: 01683659. DOI: 10 . 1016 / j . jconrel . 2011 . 10.006. URL: https://linkinghub.elsevier. com/retrieve/pii/S0168365911009588 .
• [4] Comsol. “Combining Convection and Diffusion Ef- fects”. In: (2022). URL: https : / / www . comsol . com / multiphysics / convection – diffusion – equation.

Appendice

Coating

L’applicazione del coating avviene solo se la variabile 𝚌𝚘𝚊𝚝𝚒𝚗𝚐 contiene gli indici delle sfere sulle quali applicare l’interfaccia di tipo Thin Diffusion Barrier.

if not(sum(coating==0)) surface_with_coat =[];
 for i=1:length(coating)
   surface_with_coat=[surface_with_coat , sph_boundary_alt(:,coating(i))’];
 end
 model.component(’comp1’).physics(’tds’).create(’tdb1’, ’ThinDiffusionBarrier’, 2);
 model.component(’comp1’).physics(’tds’).feature(’tdb1’).set(’ds’, num2str(d_ctg));  
 model.component(’comp1’).physics(’tds’).feature(’tdb1’).setIndex(’Ds’, D_ctg, 0);
 model.component(’comp1’).physics(’tds’).feature(’tdb1’).selection.set(surface_with_coat); % graphics to hihglight coated spheres
 figure;
 mphviewselection(model,’geom1’,surface_with_coat ,’boundary’,’facealpha’,0.5,’
 facecolorselected’,[0 1 0]) title(’Coated sphere’)
end

Integrale della concentrazione

L’estrazione della massa tramite l’integrazione della concentrazione richiedere la routine mphint2.

mass =mphint2(model,’c’,’volume’,’selection’,1);
figure; plot(initValue:step:finalValue ,mass);
title(’Mass over time’); xlabel(’Time [s]’); ylabel(’Mass [mol]’)

Linea di taglio

L’estrazione della concentrazione lungo una linea richiede la costruzione della linea e l’estrapolazione dei risultati su di essa tramite mphinterp.

for i=1:length(C_sph(1,:))
 % create cutline to read the data
 cutLine=strcat(’cln’,num2str(i));
 model.result.dataset.create(cutLine, ’CutLine3D’);
 % syntax ’genpoints ’, value , point index , coordinate
 model.result.dataset(cutLine).setIndex(’genpoints’, C_sph(1,i), 0, 0);  model.result.dataset(cutLine).setIndex(’genpoints’, C_sph(2,i), 0, 1);
 model.result.dataset(cutLine).setIndex(’genpoints’, C_sph(3,i), 0, 2);  model.result.dataset(cutLine).setIndex(’genpoints’, C_sph(1,i)+R_sph, 1, 0);
 model.result.dataset(cutLine).setIndex(’genpoints’, C_sph(2,i), 1, 1);  model.result.dataset(cutLine).setIndex(’genpoints’, C_sph(3,i), 1, 2);
 % extract results
 [conc, space,time_plot,unit]=mphinterp(model,{’c’,strcat(’cln’,num2str(i),’x’),’t’},’dataset’
,cutLine);
 [space,index]=sort(space’,’ascend’); % the data are rearranged to have the correct x distance
 x_to_plot=space; y_to_plot=conc’; % and the same index is used to  rearrange the corresponding concentration values
 for j=1:length(y_to_plot(1,:)) 
  y_to_plot(:,j)=y_to_plot(index(:,j),j);
 end
 figure; plot(x_to_plot ,y_to_plot)
 xlabel(unit{2}); ylabel(unit{1});  legend(string(time_plot(:,1)’),’Location’,’eastoutside’)  title(strcat(’Radial concentration for sphere n.’,num2str(i)))
end