Effect of chest physiotherapy on cystic fibrosis sputum nanostructure: an experimental and theoretical approach

Cystic fibrosis (CF) is a disease characterized by the production of viscous mucoid secretions in multiple organs, particularly the airways. The pathological increase of proteins, mucin and biological polymers determines their arrangement into a three-dimensional polymeric network, affecting the whole mucus and impairing the muco-ciliary clearance which promotes inflammation and bacterial infection. Thus, to improve the efficacy of the drugs usually applied in CF therapy (e.g., mucolytics, anti-inflammatory and antibiotics), an in-depth understanding of the mucus nanostructure is of utmost importance. Drug diffusivity inside a gel-like system depends on the ratio between the diffusing drug molecule radius and the mesh size of the network. Based on our previous findings, we propose the combined use of rheology and low field NMR to study the mesh size distribution of the sputum from CF patients. Specifically, we herein explore the effects of chest physiotherapy on CF sputum characteristic as evaluated by rheology, low field NMR and the drug penetration through the mucus via mathematical simulation. These data show that chest physiotherapy has beneficial effects on patients, as it favourably modifies sputum and enhances drug penetration through the respiratory mucus. Graphical abstract


Introduction
Cystic fibrosis (CF) is one of the most common lethal genetic diseases in people of Caucasian origin [1,2]. The disease is caused by mutations in the gene encoding the cystic fibrosis transmembrane conductance regulator (CFTR). The CFTR gene encodes an ATP-regulated chloride channel present within the apical surface of epithelial cells. Dis-functional CFTR determines a decreased surface liquid volume and increased mucus viscosity in many organs, mainly to the airway. In the healthy lungs, inhaled particles are transported out of the lung by cilia movement, within a process called muco-ciliary clearance. The increased viscosity of mucus in CF patients impairs muco-ciliary clearance determining mucus stagnation, thus favouring bacterial lung infections, a typical event in CF patients [2]. In this regard, CF patients are often colonized by S. aureus in the early childhood [3]. S. aureus is replaced during disease progression by P. aeruginosa as the dominant pathogen. By the time of reaching adulthood, up to 80% of CF patients are chronically infected by P. aeruginosa. These chronic pulmonary infections are associated with a rapid decline in lung functions, high morbidity and short life expectancy [4]. One of the most striking hallmarks of P. aeruginosa in chronically infected CF patients is the conversion to a mucoid phenotype of the mucus due to the increased production of alginates [5][6][7], natural polymers composed of mannuronic and guluronic acid units arranged in a block-wise pattern [8]. Alginates overproduction can contribute to a strong inflammatory response sustained by immune cells that can influence lung epithelial cells [9][10][11]. During inflammation, substances like proteins, mucin and biological polymers increase in mucus as well as the increased alginate production. These compounds organize into a three-dimensional polymeric network, affecting the whole mucus, impairing the muco-ciliary clearance, promoting bacterial infection and hindering drug diffusion to the epithelial cells.
The nanostructure of the network requires an in-depth understanding, in order to improve the efficacy of the drugs usually employed in CF therapy (e.g., mucolytics, antiinflammatory and antibiotics) when administered alone or in combination with chest physiotherapy. Indeed, drug diffusivity inside a gel-like system depends on the ratio between the diffusing drug molecule radius and the network mesh size. This parameter has driven several research studies focus to date, due to its relevance on the CF sputum/mucus properties [12]. Hanes et al. demonstrated that sputum/mucus is a heterogeneous medium made up of mainly physical crosslinked polymeric regions, embedded in a low viscosity, water-like and fluid [12][13][14]. According to their complex structure, both the macro-and the micro-rheological mucus/ sputum properties have been studied by classical strain/stress controlled rheometers and an innovative particle tracking approach [13]. The network mesh size was determined relaying on rheological sputum characteristics [14] and probe particles tracking [15]. However, as clearly pointed out by Bhat and co-workers [16][17][18], drug diffusion through mucus can be also affected by drug binding to the glycoproteins. In addition, drug action can be reduced by non-physiological ionic strength, pH, mucociliary clearance and biochemical alterations of mucus [18].
Based on our previous findings, we herein propose the synergic combination of rheology and low-field nuclear magnetic resonance (LF-NMR) to study the mesh size distribution changes of mucus in CF patients (obtained by voluntary expectoration and from now on called "sputum"), prior and after chest physiotherapy. Indeed, while LF-NMR can provide for info about sputum structure, depending on the magnetic relaxation of water hydrogen trapped in the network, rheology can yield to further insights on the nanostructure, relaying on the mechanical relaxation of the polymeric network chains. Furthermore, LF-NMR allows determining the spin-spin relaxation time (T 2m ) of the water hydrogens, whereas rheology allows measuring the shear modulus G of the sputum. Remarkably, it is possible to determine important features of the sputum nanostructure and further develop innovative drug penetration range, on the T 2m and G knowledge. In particular, the present work focuses on the effects that chest physiotherapy has on sputum characteristic leant on rheology/LF-NMR analysis. Therefore, we hypothesize that chest physiotherapy can also influence the sputum nanostructure, representing a pivotal information for clinicians in order to optimize the therapy, for future precision and personalized medicine approach.

Samples collection
Sputum samples were provided by the Burlo Garofolo Hospital, following a procedure approved by the Ethics Committee (prot n. 496/2916, CI M-11, 22-3-2016). Written informed consent was obtained from each patient. One sputum specimen was collected before and after chest physiotherapy session in randomly selected adult and young patients able to expectorate. Non-expectorant patients were excluded, as were subjects ≤ 7 years of age. The chest physiotherapy technique is set up on patient to patient even if the technique most used for the patients enrolled in this study is the pep-mask. Sometimes, autogenous drainage has also been used. Chest physiotherapy was provided as standard procedure for CF patients, to improve clearance, both by using positive expiratory pressure devices and standardised physiotherapist-guided chest applied manoeuvres, as well as effective breathing exercise [19].
The aforementioned criteria allowed for the recruitment of most patients who attended Burlo Garofolo Hospital for a clinical visit during our study period (from January 2019 to December 2020). Spontaneously expectorated (1-2 mL) sputum was collected from CF patients in sterile cups and immediately used for T 2m determination. The samples were then transferred from the LF-NMR glass tube to a rheometer device for analysis. Patient's lung functionality, performed after chest physiotherapy, was evaluated by the Burlo Garofolo Hospital by means of FEV 1 (forced expiratory volume in the first second), i.e. the volume of air that can be forced out in one second after taking a deep breath.

Rheology
The viscoelastic properties of the analysed samples were determined by means of two rheological tests called, respectively, stress sweep and frequency sweep test. These experiments require the application of a sinusoidal stress ( ) to the sample in order to measure the related oscillatory deformation ( ). Stress sweep tests were lead applying a sinusoidal stress of constant frequency (f = 1 Hz) and increasing amplitude to detect the limit of the linear viscoelastic range. Frequency sweep tests, performed in the linear viscoelastic range, as to prevent any sample structure damaging, were lead in the frequency range 10-0.01 Hz. The output of the frequency sweep tests, i.e. the elastic (G') and viscous (G'') moduli dependence on pulsation (= 2 f) (mechanical spectrum), were fitted by the generalized Maxwell model [20]: where n R is the number of Maxwell elements considered, whereas g i (ith elastic constant), i (ith viscosity) and i (ith relaxation time) represent model fitting parameters. Model fitting was performed assuming that i were scaled by a factor 10 ( i+1 = 10 i ) [20]. Although other fitting strategies could have been followed, this approach enables getting a wider and, thus, clearer, relaxation spectrum (g i vs i ). The n R determination was performed according to a statistical procedure, in order to minimize the product N* 2 , where N indicates the number of fitting parameters and 2 is the sum of the squared errors [21]). The sample shear modulus G was evaluated as the sum of all g i (G = Σg i ) [22]. Therefore, Flory theory [23] allows to evaluate the polymeric network crosslink density, x , defined as the moles of junction points between different polymeric chains per hydrogel unit volume, based on previous G determination: where R is the universal gas constant and T is the absolute temperature. Therefore, the equivalent network theory [24] enables the determination of the average mesh size a : where N A is the Avogadro number.

Low-field NMR
NMR is based on hydrogen atoms dipole reaction to the external constant magnetic field (B 0 ) switch where they are embedded. Indeed, inside B 0 , hydrogen atoms dipoles tend to line up in the B 0 direction so that, globally, they give origin where t is time, I(t) is the ratio between the time dependent value of M XY and its maximum value (M XYmax ) occurring just after B 1 removal, T 2i represent the ith spin-spin or transverse relaxation times, while A i are "weights" proportional to the number of hydrogen atoms whose dipoles relaxation is characterized by T 2i . Equation (5) simply states that the relaxation process is the result of "m" exponential relaxation processes each one characterized by its own relaxation time (T 2i ) and weight (A i ). The determination of the unknown couples (T 2i , A i ) was performed by Eq. (5) fitting to the experimental I(t) values (I s (t)), and the number, m, of exponential decays appearing in Eq. (5) was determined per the statistics applied in the Rheology section [21]. The average relaxation time (T 2m ) and the average inverse relaxation time ((1/T 2 ) m ), depending on several variables such as temperature, B 0 strength, and the presence of a disperse phase in the system as it happens in hydrogel system [26], can be defined by: when T 2m is evaluated as the inverse of (1/T 2 ) m , it is named T 2mb .
While Eq. (5) provides the discrete relaxation time distribution represented by the m couples (A i% -T 2i ), it is possible to determine the continuous distribution according to what suggested by Whittal and MacKay [27]: where T 2max (= 10 4 ms) and T 2min (= 0.1 ms) indicate, respectively, the lower and upper values of the continuous T 2 distribution, a(T 2 ) is the unknown amplitude of the spectral component at the relaxation time T 2 , while exp{-t/T 2 } represents the decay term. Equation (7) represents the "continuous" expression of Eq. (5). This means that while Eq. (5) describes the M XY relaxation process by means of a limited (discrete) number of initially unknowns T 2i , Eq. (7) makes use of a much higher number of known relaxation times belonging to the interval (T 2min -T 2max ) and comprehending all the possible relaxation times characterizing the sample under study (continuous distribution of relaxation times). In order to fit the experimental M XY time decay (I s (t)) by Eq. (7) and get the continuous T 2 distribution (the unknowns A i = a i (T 2i )*ΔT 2i ), the following discretization was applied [27]: where the range of the T 2 distribution (T 2min -T 2max ) was logarithmically subdivided into N = 200 parts (higher N values were unnecessary). Because of the noise disturbing the measure of I s , the fitting procedure must not minimize the 2 statistic, but a smoothed definition of it ( 2 s ) [27]: where i is the i th datum standard deviation and is the weight of the smoothing term (second summation in Eq. (9)) proposed by Provencher [28]. Although different criteria can be used to determine , the strategy proposed by Wang [29] was applied. Based on this strategy, the correct value is that occurring just below the heel (slope variation) of the function ln( s ) vs ln( ). In the present work, = 150 was determined.
The discrete and continuous T 2 distribution can be transformed into hydrogel mesh size distribution resorting to one of the fundamental relations of the low-field NMR field. This relation, based on the solution of the magnetization diffusion equation proposed by Brownstein and Tarr [25,26], establishes the link between (1/T 2 ) m and the ratio of the surface (S) of the dispersed/solubilized substances in the sample and the volume (V) of the sample water molecules: where T 2H2O is the bulk protons relaxation time (i.e. the water proton relaxation time in the absence of polymer, the so-called free water relaxation time ≈ 3700 ms at 37 °C and B 0 = 0.47 T [30]) and M (length/time) is a physical parameter, named relaxivity, which represents the effect of the surface of polymer chains on water proton relaxation. Indeed, M is equal to the ratio between thickness and relaxation time of the bound water layer adhering to the solid surface. Equation (10), stating that (1/T 2 ) m depends on (S/V), clearly establishes the relation between the relaxation time and the spatial organization of the sample network that heavily affects the S/V ratio [26]. For instance, in many polymeric systems, crosslinking induces a spatial reorganization of the polymeric chains contained in the original solution that involves the increase of the ratio S/V [31,32]. This, in turn, reflects in the increase of (1/T 2 ) m and in the decrease of T 2m . Despite its theoretical importance, Eq. (10) can be rewritten in a more useful form, based on the fibre cell [26] and Scherer theories [33]. At this purpose, Abrami et al. [34] demonstrated that, for a hydrogel polymer volume fraction < 0.6, the term (S/V) can be expressed as a function of the average mesh size a of the polymeric network (error < 5%): where C 0 and C 1 are two constants depending on the mesh architecture and equal, respectively, to 1 and 3 in the case of cubic mesh [33]. Thus, Eq. (10) becomes while Eq. (12) refers, averagely, to all the polymeric network meshes, similar expressions can be written for meshes of different dimensions ( i ), assuming the M independence on the mesh size [26]: where T 2i is the relaxation time of water protons trapped in polymer meshes of size i . The bi-univocal correspondence between T 2i and i only holds in the fast-diffusion regime (typical of gels), i.e. when the mobility of water molecules, expressed by their self-diffusion coefficient D (3.04 × 10 −9 m 2 /s at 37 °C [35]), is large whether compared to the rate of magnetization loss (R c × M) (i.e. R c *M/D ≪ 1). In the slow diffusion regime, relaxation of all the water protons contained in the volume of a mesh of size i is not described by just one T 2i but several T 2i . R c indicates the radial distance from the polymer chain axis where the effect of polymeric chains on water proton relaxation becomes negligible. This can be expressed by [26] The combination of Eqs. (12) and (13) allows to conclude that the ratio between i and its average value, a , depends exclusively on the relaxation times T 2i and T 2m (except for the free water relaxation time T 2H2O ): Thus, Eq. (5) or Eq. (8) fitting to experimental M XY decay (I s (t)) allows determining the discrete or continuous relaxation spectrum (A i% − T 2i ), whereas Eq. (13) provides for the conversion into the discrete or continuous mesh size distribution of the sputum (A i% − i ).
LF-NMR measurements were performed at 37 °C by means of a Bruker Minispec mq20 (0.47 T, Germany). The determination of T 2m was performed according to the CPMG sequence (Carr-Purcell-Meiboom-Gill) [36] {90°[-τ-180°τ(echo)]n-T R } with a 8.36 μs wide 90° pulse, τ = 250 μs, and T R (sequences repetition rate) equal to 10 s. In order to obtain the final I s of 2%, the proper n was chosen. m was determined according to the statistics applied in the Rheology section [21]. Each spin-echo decay, composed by n points, was repeated 36 times (number of scans).

Drug diffusion
The walls of the bronchial tree are covered by a very thin layer of a water-like serous fluid, the periciliary liquid (PCL), where cells cilia beat with a typical frequency of 20 Hz and amplitude about 5 m. On the top of PCL, a thin layer of a non-Newtonian viscoelastic fluid lies, named mucus [37] (see Fig. 1).
PCL and mucus layer (ML) constitute the airway surface liquid (ASL) [38]. The viscoelastic nature of mucus depends on the presence of different components such as proteins, alginates, white blood cells, DNA, bacteria and mucin (whose levels are pathologically increased compared with healthy subjects [39]) that lead to a three-dimensional network [40]. Drug permeation inside the ASL is a complex phenomenon affected by many factors such as the thickness of ASL, which depends on the ASL rheological properties and the airflow [37,41,42]. Moreover, possible interactions (e.g. electrostatic and hydrogen bonding) between the diffusing drug and the polymeric components of the ALS cannot be a priori excluded [43]. Furthermore, as ASL is not homogeneous, the drug diffusion coefficient (D) should be position dependent [42]. However, assuming that drug transport is mainly affected by diffusion with constant D, the prediction of the drug penetration inside ASL can be obtained by referring to the one-dimensional Fick's equation with proper initial and boundary conditions. In detail, it is assumed that at time zero (deposition of drug solution on ASL, following inhalation), the drug is not present inside ASL, and it is uniformly distributed, at the known concentration C 0 , in the solution layer of thickness (see Fig. 1). Boundary conditions involve the existence of an impermeable wall in X = 0 (the drug cannot diffuse in the X negative direction) and the possibility of drug diffusion in the positive X direction up to infinite (therefore, the drug can diffuse inside cells, too). According to these hypotheses, the time evolution of the drug profile concentration (C) inside ASL is determined by where t is time. The D value to be used in Eq. (15) can be provided by the Lusting-Peppas equation [44] as per Abrami et al. [23]: where r s is the radius of the drug molecule imagined to be a sphere, D 0 is the drug diffusion coefficient in water, while Y is a model parameter that, lacking further information, can be set equal to 1 [44], although it usually ranges between 3 and 30 [45,46] (see Appendix 1 for some considerations on the potentiality of Eq. (17)). The solid volume fraction φ was set equal to 0.05, this being the most common situation, albeit variability can occur among samples [47]. For what concerns r s and D 0 , with an attention focus on typical antibiotics used for FC patients (ciprofloxacin and tobramycin, for example [48]), reasonable values are 1.1 nm and 3 × 10 −10 m 2 /s, respectively [49,50]. Moreover, it is worth mentioning that while the PCL thickness ranges between 5 and 10 m [51], in healthy subjects and in CF patients, ML thickness is similar and, approximately, equal to 25-30 m [37,52]. Vice-versa, patients affected by chronic obstructive pulmonary diseases (COPD) exhibit an ML thickness value up to 300 m [41].

Statistical analysis
The nature of the experimental data distribution (normal or not) was evaluated by the Kolmogorov-Smirnov test (KStest). Based on KS-test results, the Spearman's correlation coefficient (r sp ) was considered to verify possible direct or inverse correlations among T 2m , FEV 1 and G. Correlations among variables were considered significant when p <0.05, corresponding to a probability of 95%. Lower probability was associated to a lack of correlation among variables. Figure 2 remarkably shows that there is a statistically significant correlation between T 2m measured in CF sputum and FEV 1 before (dark circles; r sp = 0.69, p = 0.006) and after (white circles; r sp = 0.44, p = 0.044) chest physiotherapy. As lung functionality increases with FEV 1 [53], Fig. 2 indicates that T 2m increase is connected with improved patient clinical conditions. Moreover, the inspection of Fig. 2 reveals that, on average, after chest physiotherapy, T 2m is higher. Indeed, the dashed line (connected with the T 2m vs FEV 1 trend after chest physiotherapy) is characterized by a higher intercept with respect to solid line (connected with the T 2m vs FEV 1 trend after chest physiotherapy), and the two straight lines share, approximately, the same slope. Thus, the beneficial action of chest physiotherapy is clearly evident.

Results and discussion
Interestingly, similar correlations between T 2mb (average relaxation time evaluated as the inverse of (1/T 2 ) m ) and FEV 1 have been determined both before and after chest physiotherapy as shown on Fig. 3. This finding is due to the strong correlation existing between T 2m and T 2mb (r sp = 0.97, p <0.0001), as per Fig. 4  for all the analysed samples. Even though considering separately the samples before and after chest physiotherapy, the same strong correlation is demonstrated too.
Thus, for what concerns the evaluation of lung functionality in relation to FEV 1 , both T 2m and T 2mb can be chosen as indicators of patient clinical conditions. Figure 5 shows the significant inverse correlation between T 2m and the shear modulus G on the samples both before and after chest physiotherapy (r sp = −0.58, p = 0.0001). These evidences support the beneficial effects of chest physiotherapy. Indeed, at fixed G, Fig. 5 shows that chest physiotherapy leads to a higher T 2m (see the solid and dashed lines in Fig. 5). An additional and conclusive proof of chest physiotherapy benefit is displayed in Fig. 6. This picture reports the %, relative (compared to the value before chest physiotherapy, T 2mBef ) variation ( ΔT 2m ) of T 2m versus the %, relative (with respect to the value before chest physiotherapy, G Bef ) variation ( Δ G) of G. It is possible to detect that in the 57.1% of the cases (section I in Fig. 6), chest physiotherapy leads to both an increase of T 2m and a decrease of G, which indicates the partial breakage of the original mucus polymeric network that yields to a weaker nanostructure. Indeed, the combination of Eqs. (3), (4) and (12) evidences that T 2m (or T 2mb ) is inversely proportional to G 1/3 . Therefore, network damaging is due to a T 2m increase, as per G decrease due to a reduction of crosslink density (see Eq. (3)). The direct correlation between T 2m and FEV 1 evidences that nanostructure variation indicates a better lung functionality. Figure 6, section IV shows that in the 23.9% of the cases, chest physiotherapy determines an increase both in T 2m and G. This behaviour could be explained assuming that the increase of sputum stiffness is not due to the increase of the crosslink density, but to the increase of polymeric fibres strength as the polymeric chains group into thicker fibres (see Fig. 6). This perfectly matches with G and T 2m increase, as T 2m rises with the decrease of the ratio between the solid (fibres) surface (S) in contact with water molecules and the hydrogel water molecule volume (V) as per Eq. (10) [22,54].  Fig. 6). Depending on the relative variation of T 2m and G, different mucus nanostructures can be supposed Figure 6 section II reveals that in 9.5% of the cases, chest physiotherapy determines both G and T 2m reduction. This proves the increase of the crosslink density, which takes place among thinner and weaker fibres. In the 9.5% of the cases, chest physiotherapy determines a G increase and T 2m decrease ( Fig. 6 section III). In this case, T 2m decrease and G increase could be related to mesh size reduction. Overall, the effect of chest physiotherapy on sputum nanostructure is positive, as it determines the average T 2m increase and the average G decrease (see the black dot in section I of Fig. 6). Figure 7 qualitatively sums up the effect of chest physiotherapy on the sputum nanostructure, based on the relative variations of T 2m and G shown in Fig. 6.
It is of outmost importance to remind that all the above considerations relay on the interpretation of experimental results in the light of Eqs. (3), (4) and (12). Thus, other experimental techniques (such as TEM/SEM or other image techniques) will be considered. However, the complex nature of CF sputum nanostructure, reported in Fig. 8, and the issues connected to the reliable acquisition and interpretation of images [55][56][57] make this task absolutely challenging.
Moreover, for what concerns the image analysis, the presence of DNA, cell debris and bacteria can create an extremely complex scenario. On the other hand, this is not a problem for the proposed combined approach (rheology and LF-NMR) as the presence of additional components just alters the S/V ratio (see Eq. (10)) that reflects in a variation of T 2m . Similarly, the presence of additional components affects the sputum rheological behaviour reflecting in alterations of G' (shall they represent physical or chemical bonds among different polymeric chains) or G'' (shall they stay in the sputum sol fraction). Thus, the presence of DNA, cell debris and bacteria is accounted for by our approach, and it is interpreted as a sort of network architecture modifications.
The general soundness of the presented approach is proved by Figs. 9 and 10 showing, respectively, the average mesh size ( a ) of all samples and the continuous mesh size distribution referring to sample 11 before and after chest physiotherapy (see Table 1 in Appendix 2). Indeed, according to Hanes and co-workers [58,59], who applied a sophisticated particles tracking methods, the mesh size distribution of the CF sputum network ranges between 60 and 300 nm, with a mean value of (145 ± 50) nm. Figure 9 shows that, except for one case, a refers to the range of 50-250 nm, i.e. within the range indicated by Hanes and co-workers. Moreover, Fig. 10 shows that before chest physiotherapy (unperturbed sputum), the continuous mesh size distribution spans between 20 and 300 nm, perfectly matching Hanes and co-workers findings. Thus, the two different strategies (our bulk rheology-LF-NMR and Hanes' particle tracking) provide similar results. Figure 9 also allows to sum up the effect of chest physiotherapy on sputum nanostructure. Indeed, Fig. 9 suggests that, on average, chest physiotherapy leads to an increase of a , as its relative variation (100Δ a / bef ) is mainly positive. Conversely, negative variations are limited in their amplitude (grey circles), with reference to positive ones (see dotted grey line). As per Eqs. (8) and (15), it is possible to determine the variation of the continuous mesh size distribution,  Figure 10 shows, e.g. this comparison in the case of samples 11 (see Table 1 in Appendix 2).
It is possible to notice that chest physiotherapy pushes towards higher mesh size values ( ) the entire mesh size distribution (compare the black solid line and dashed black line). It also provides for changes in the shape of the distribution, making more evident the second peak, which is just slightly visible before the treatment. This is a pivotal hallmark: the second peak, moving from, about, 200 nm to 500 nm, refers to mesh size range which is typical of healthy sputum.
Moreover, Fig. 10 clearly proves that about 20% of the mesh size distribution is on the normal range size (dashed black line), while, before chest physiotherapy (continuous solid line), 100% of the mesh size refers to the pathological range. Therefore, both before and after chest physiotherapy, the discrete mesh size distribution (grey vertical lines) results in a trend of the continuous one.
Drug delivery through mucus is a complex phenomenon which is on stage for researchers, with a main focus of the current COVID period the world is facing [16][17][18]60]. For what concerns drug (e.g. antibiotics, anti-inflammatory and mucolytics) administration by inhalation and in order to achieve this task, detailed information about mesh size are of key importance. Indeed, the entire delivery process implies different steps such as drug deposition on ASL (airway surface liquid, sum of the periciliary liquid and mucus layers) and drug diffusion throughout it. Modern strategies ensure that within, about, 7 min, the administered aerosol can get the ASL [61,62]. In the light of the simplifying hypotheses on which Eq. (16) relies, we would like to simulate the kinetics of the drug diffusion process inside ASL. At this purpose, the drug diffusion coefficient D was evaluated according to Eq. (17) assuming the values of the model parameters (φ = 0.05, r s = 1.1 nm, D 0 = 3 × 10 −10 m 2 /s) presented in the "Diffusion" section. Moreover, the average mesh size was fixed equal to the average value (92 nm) reported in Table 1; Y was set equal to 3 or 30 to account for two limiting conditions and = 3 m. When Y = 3 (this corresponding to a higher drug diffusivity D = 2.5 × 10 −10 m 2 /s), Fig. 11A reveals that within few seconds, the drug can cross the ASL typical thickness of healthy subjects and CF patients (≈ 30 m).
Indeed, the dimensionless concentration C/C 0 achieves about 50% of the initial solution concentration at the ASL bottom in a very short time (see the profile concentration curves labelled 1.2 and 3 s in Fig. 11A). For longer times, the C/C 0 profile assumes a flatter shape (see the 15 s labelled curve in Fig. 11A) and a uniform distribution is achieved after 240 s through the ASL thickness competing to COPD (chronic obstructive pulmonary diseases) patients (300 m). COPD has been considered as it is due to an acquired  Fig. 10 Comparison between the continuous mesh size distribution before (solid black curve) and after (dashed black line) chest physiotherapy referring to samples 11 (see Table 1). Solid grey vertical lines (secondary vertical axis on the right) indicate the discrete mesh size distribution referring the solid black curve. Dashed grey vertical lines (secondary vertical axis on the right) indicate the discrete mesh size distribution referring the dashed black curve CFTR impairment, resulting in increased lung mucus viscosity similarly to CF. Considering Y = 30, this value represents a lower value of the diffusion coefficient (D = 0.6 × 10 −10 m 2 /s). Figure 11B reveals that the diffusion process is slowed down, but C/C 0 is about 50% at a penetration depth of 30 m after 5 s. Moreover, 900 s are required to get a flat concentration profile up to 300 m, i.e. the ASL thickness which can be detected in COPD patients. Thus, in the case of CF patients, these simulations indicate that the drug can diffuse on the ASL thickness in a few seconds. On the other hand, in the case of COPD patients, 6-15 min are required to get the same effect if ASL thickness is about 300 m. It is worth to mention that the values of the diffusion coefficients therein reported, assuming Y = 3 or 30 (2.5 × 10 −10 m 2 /s and 0.6 × 10 −10 m 2 /s, respectively), approximately represent the maximum and minimum D range which can be determined by the experimental permeability data through mucus performed by Bhat and co-workers [17,18] that span from 2.4 × 10 −10 m 2 /s down to 0.7 × 10 −10 m 2 /s. In the case of CF patients, drug diffusion through mucus is not the rate determining step of the whole delivery process. In the worst conditions for COPD patients, drug diffusion can be compared to the inhalation time. However, in both cases, the time required for drug diffusion is short, this witnessing the promising role of inhalation to match pulmonary diseases.

Conclusions
The characterization of the sputum of CF by the combined use of rheology and LF-NMR allowed getting information about the nanostructure of sputum derived from the mucus lining the bronchial tree. In particular, this analysis was aimed to determine sputum nanostructure variations, as a result of chest physiotherapy, a routine procedure adopted in CF patients. We could verify that significant nanostructure variations occurred upon chest physiotherapy, leading to beneficial clinical effects for the patients. Interestingly, the combined use of rheology and LF-NMR allowed to define the nanostructure characteristics that each technique (if used singularly) would have never been able to provide.
Moreover, assuming that drugs diffuse inside mucus according to Fick's law, as well as in the sputum nanostructure, we could determine that in the case of healthy subjects and CF patients, about 7 min from inhalation are needed to get the drug active inside the mucus, i.e. a very short time. In case of COPD patients, this time has been estimated to about 15 min in the worst situation (300 m ASL thickness). Anyway, in both cases, the time required to get drug action since inhalation resulted to be short. Therefore, the inhalation delivery demonstrated a key role to provide for effective info about pulmonary diseases, such as those afflicting CF and COPD patients. Funding Open access funding provided by Università degli Studi di Trieste within the CRUI-CARE Agreement. This paper has been financed by Fondazione CRT Trieste project n°103069 and by the socalled Programma di valorizzazione dei brevetti del sistema universitario del Friuli Venezia Giulia, Unity FVG PoC, 2020, Italy.

Availability of data and materials
Where not directly provided in the paper, all the data referring to this paper are available upon request to the corresponding author.

Declarations
Ethics approval and consent to participate Sputum samples were provided by the Burlo Garofolo Hospital, following a procedure approved by the Ethics Committee (prot n. 496/2916, CI M-11, 22-3-2016). Written informed consent was obtained from each patient.

Competing interests
The authors declare no competing interests.

Appendix 1
The approximation of the Scherer theory [33] by Abrami and co-workers [34] allows to provide a simple relation between the average mesh size a , the polymer volume fraction φ and the polymeric chains radius r f : where C 1 and C 0 are two constants depending on the network architecture (cubical, tetrahedral or octahedral) whose values are defined by the Scherer theory ( [33]; for a cubical network, C 0 = 1 and C 1 = 3 ). The combination of Eq. (18) and Eq. (17) (19) assuming different values of the ratio between the drug molecule radius (r s ) and the radius (r f ) of the fibres constituting the polymeric network pervading the hydrogel. These simulation refer to a cubical mesh organization (C 0 = 1; C 1 = 3 [33])  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.