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

Ovvero 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 Renolds.

\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 inter- agire 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.

FIG. 2: Render grafico per la descrizione del problema affrontato. Sono presenti quattro sfere cariche di farmaco all’interno di un cilindro contenente una soluzione liquida e si vuole analizzare la cinetica diffusiva che si instaura tra le sfere, inizialmente cariche, e la soluzione.
FIG. 2: Render grafico per la descrizione del problema affrontato. Sono presenti quattro sfere cariche di farmaco all’interno di un cilindro contenente una soluzione liquida e si vuole analizzare la cinetica diffusiva che si instaura tra le sfere, inizialmente cariche, e la soluzione.

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.

VariabileValoreSignificato
R_cyl3E-2Raggio del cilindro
H_cyl2E-2Altezza del cilindro
R_sph0.3E-2Raggio delle sfere
D_liq1E-6Coefficiente di diffusione per la soluzione
D_sph2E-6Coefficiente di diffusione per le sfere
C03Concentrazione iniziale
Y_ctg0.1E-3Ammettenza del coating
d_ctgR_sph/1000Spessore del coating
TAB. 1: Parametri materiali usati nella formulazione del problema

Il codice parte caricando i pacchetti Comsol e attivando il percorso verso la cartella conte- nente le subroutines. Dopo la prima fase di in- izializzazione 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 ri- aprire 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 mod- ello 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) precedente- mente 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.

FIG. 3: Diagramma a blocchi descrittivo della rou- tine principale per l’impostazione e la risoluzione di tale problema diffusivo. Il loop tratteggiato 𝚖𝚎𝚜𝚑 𝚜𝚎𝚗𝚜𝚒𝚝𝚒𝚟𝚒𝚝𝚢 viene sfruttato solo inizialmente nell’analisi di sensibilità al magliaggio. Una volta selezionato il grado di mesh viene completata anche l’analisi considerando il coating.
FIG. 3: Diagramma a blocchi descrittivo della routine principale per l’impostazione e la risoluzione di tale problema diffusivo. Il loop tratteggiato 𝚖𝚎𝚜𝚑 𝚜𝚎𝚗𝚜𝚒𝚝𝚒𝚟𝚒𝚝𝚢 viene sfruttato solo inizialmente nell’analisi di sensibilità al magliaggio. Una volta selezionato il grado di mesh viene completata anche l’analisi considerando il coating.

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 nu- merico 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 tempo- rale 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. Successiva- mente vengono estratti i dati esattamente sulle cutline e sono riportati in fig. 6b.

FIG. 4: Concentrazione lungo il raggio della sfera n. 4 a partire dal centro verso il bordo. Le curve, dall’alto, variano al variare del tempo, da 0 a 4s.
FIG. 4: Concentrazione lungo il raggio della sfera n. 4 a partire dal centro verso il bordo. Le curve, dall’alto, variano al variare del tempo, da 0 a 4s.

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.

FIG. 5: Risultati dell’analisi di convergenza (a); mesh a livello Extra Fine (b). Sono raffigurati in scala loga- ritmica i gradi di libertà. Sulle ordinate sono presenti l’errore percentuale (in blu) e il tempo di simulazione (in rosso). Appare evidente come l’errore relativo si riduce sotto il 5% con la mesh di ID 2 dove il tempo è 5 volte inferiore alla discretizzazione massima.
FIG. 5: Risultati dell’analisi di convergenza (a); mesh a livello Extra Fine (b). Sono raffigurati in scala loga- ritmica i gradi di libertà. Sulle ordinate sono presenti l’errore percentuale (in blu) e il tempo di simulazione (in rosso). Appare evidente come l’errore relativo si riduce sotto il 5% con la mesh di ID 2 dove il tempo è 5 volte inferiore alla discretizzazione massima.

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).

Mesh refinementExtr. coarseExtra coarseCoarserCoarseNormalFineFinerExtra fineExtr. Fine
ID987654321
DoF797252044938924130441969745692154456718178
Error8.22201.15910.51910.22230.22150.20360.09130.02140
Time cost5.98204.52294.41235.09825.83786.436810.113827.1032133.6901
TAB. 2: Risultati dell’analisi di convergenza. Viene riportato l’ID numerico rappresentante in Comsol la mesh selezionata, il numero di gradi di libertà, l’errore normalizzato e il costo computazionale in secondi.

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).

FIG. 6: Variazione della massa nel tempo nel cilindro (a) e nelle sfere (b). Tali curve rappresentano la quantità di
farmaco che diffonde al passare del tempo.
FIG. 6: Variazione della massa nel tempo nel cilindro (a) e nelle sfere (b). Tali curve rappresentano la quantità di farmaco che diffonde al passare del tempo.

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.

FIG. 7: Plot di superficie con la concentrazione all’instante iniziale (a) e finale (b)
FIG. 7: Plot di superficie con la concentrazione all’instante iniziale (a) e finale (b)

Coating superficiale

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

FIG. 8: Geometria di riferimento con numerazione delle sfere (a) dove 1,2,3,4 corrispondono agli ID Comsol di dominio 2,3,4,5. Sfere sulle quali è stato applicato il coating (b)
FIG. 8: Geometria di riferimento con numerazione delle sfere (a) dove 1,2,3,4 corrispondono agli ID Comsol di dominio 2,3,4,5. Sfere sulle quali è stato applicato il coating (b)

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).

FIG. 9: Variazione della massa nel tempo nel cilindro (a) e nelle sfere (b). Si vede nettamente la differenza tra le sfere con coating (dominio 2,5) e senza coating (3,4).
FIG. 9: Variazione della massa nel tempo nel cilindro (a) e nelle sfere (b). Si vede nettamente la differenza tra le sfere con coating (dominio 2,5) e senza coating (3,4).

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.

FIG. 10: Andamento radiale della concentrazione nelle sfere con coating (a) e senza coating (b) al variare del tempo di simulazione. La coordinata orizzontale si muove dal centro della sfera al bordo.
FIG. 10: Andamento radiale della concentrazione nelle sfere con coating (a) e senza coating (b) al variare del tempo di simulazione. La coordinata orizzontale si muove dal centro della sfera al bordo.

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 prob- lema 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).

FIG. 11: Variazione della massa nel tempo (estensione del tempo di simulazione a 40s) nel cilindro (a) e nelle sfere (b). Risulta evidente la differenza tra le sfere con coating (domini 2,5) e senza (domini 3,4). La massa residua tende ad diventare nulla e diffonde nel cilindro ma lo fa con velocità diverse tra le due famiglie di sfere.
FIG. 11: Variazione della massa nel tempo (estensione del tempo di simulazione a 40s) nel cilindro (a) e nelle sfere (b). Risulta evidente la differenza tra le sfere con coating (domini 2,5) e senza (domini 3,4). La massa residua tende ad diventare nulla e diffonde nel cilindro ma lo fa con velocità diverse tra le due famiglie di sfere.

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 ge- ometrie 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