Olivier Lafforgue
email: olivier.lafforgue@univ-amu.fr

Isabelle Seyssiecq
email: isabelle.seyssiecq@univ-amu.fr

Sébastien Poncet
email: sebastien.poncet@usherbrooke.ca

Julien Favier
email: julien.favier@univ-amu.fr

Laurent Morin
email: morin@physio-assist.com

A time-dependent Herschel-Bulkley model for synthetic bronchial mucus

Keywords: Synthetic mucus, rheology, Herschel-Bulkley model, thixotropy

Airway mucus is a complex material from a rheological point of view exhibiting a yield stress, shear-thinning, viscoplastic and thixotropic properties. Collecting real mucus being quite challenging and inducing a huge variability of the samples, mucus simulants are usually preferred to investigate in details its rheology and derive behavioral law, which are of prime importance to develop helping clearance devices and advanced numerical models. From new measurements through series of stepwise changes of shear rate, a time-dependent Herschel-Bulkley model is here proposed for a synthetic mucus representative of cystic fibrosis patients. The analysis is decomposed into a simplified model for steady state conditions and a full unsteady model.

I. INTRODUCTION

Bronchial mucus plays a crucial role in the protection of the respiratory system and can be seen as a key-element for better understanding and treating chronic respiratory diseases such as severe asthma or cystic fibrosis. The airway clearance function is a complex phenomenon directly linked to the thermophysical and mostly rheological properties of the mucus layer coating the airways. The use of real mucus implies strong issues related to available quantities, and rises questions about the impact of the collection method on the fluid composition. Moreover, disease conditions can critically drive the mucus away from its normal physiological state and thus from its original rheological properties by affecting either mucins size, type and concentration [START_REF] Thornton | From mucins to mucus: toward a more coherent understanding of this essential barrier[END_REF]. Its huge variability led us to choose working with mucus simulants.

Former works focusing on the mucus rheology, either synthetic [START_REF] Hamed | Synthetic tracheal mucus with native rheological and surface tension properties[END_REF] or native [START_REF] Bastholm | The viscoelastic properties of the cervical mucus plug[END_REF], have only described part of its rheological properties, using dynamic oscillatory shear measurements. Under small deformations these measurements mostly reflect the properties of mucus under its unperturbed Contract grant sponsor: Physio-Assist and the Association Nationale de la Recherche et de la Technologie (grant number 20141287). state and are useful to describe the linear viscoelastic response of mucus. On the contrary, some authors [START_REF] Zahm | Role of simulated repetitive coughing in mucus clearance[END_REF], [START_REF] Lee | Effects of chelating agents on the rheological property of cervical mucus[END_REF] made the choice to characterize mucus rheology using only continuous shear experiments, which mostly reflect the flow behavior of mucus and can be used to investigate its viscosity under physiological shearing rates prevailing in human lungs during normal functioning or during temporary events such as cough. Lafforgue et al. [START_REF] Lafforgue | Thermo-physical properties of synthetic mucus for the study of airway clearance[END_REF] provided the widest experimental database of the temperature-dependent physical properties, namely surface tension, heat capacity, density, thermal conductivity and diffusivity, of the present synthetic mucus. It has been later extended by [START_REF] Lafforgue | Rheological properties of synthetic mucus for airway clearance[END_REF] to its main rheological characteristics. They showed that by changing the concentration of Actigum within the solution, they were able to mimic real mucin concentration for different respiratory pathologies, from asthma to cystic fibrosis. From a rheological point of view, this synthetic mucus exhibited a viscoplastic and shear-thinning behavior with a yield stress and hysteresis loops in agreement with the scientific literature. The empirical parameters of the underlying Herschel-Bulkley law are dependent on the Actigum's concentration. This database enabled to improve the design of a helping clearance device tested clinically in [START_REF] Giovannetti | Étude SIMETOL: évaluation de la tolérance et l'efficacité d'une procédure innovante de drainage bronchique assistée par Simeox dans la mucoviscidose en alternative à la kinésithérapie respiratoire conventionnelle[END_REF].

The thixotropic properties of airway mucus have not been considered so far. One can cite the works of Nielsen et al. [START_REF] Nielsen | Elastic contributions dominate the viscoelastic properties of sputum from cystic fibrosis patients[END_REF] and Tomaiuolo et al. [START_REF] Tomaiuolo | A new method to improve the clinical evaluation of cystic fibrosis patients by mucus viscoelastic properties[END_REF], who observed an hysteresis loop on the continuous flow curve for samples of human sputum collected from patients suffering from cystic fibrosis. The interested reader can refer to the paper of Lai et al. [START_REF] Lai | Micro-and macrorheology of mucus[END_REF], who proposed the most complete review on the mucus rheology til 2009.

The present paper proposes to extend the rheological analysis of Lafforgue et al. [START_REF] Lafforgue | Rheological properties of synthetic mucus for airway clearance[END_REF] to quantify the thixotropy of synthetic mucus and model it by a simplified then full timedependent Herschel-Bulkley model.

II. EXPERIMENTAL METHODOLOGY A. Preparation of synthetic mucus

Mucus simulants were prepared within glass bottles filled with 200 mg of distilled water. Each component was then very slowly poured into the solution stirred using a magnetic stirrer at room temperature (around 21 • C). The protocol proposed by Zahm et al. [START_REF] Zahm | Role of simulated repetitive coughing in mucus clearance[END_REF] consists in the addition, in the following order, of 0.9wt% of NaCl, 0.5wt% of Viscogum FA (galactomannan gum derived from locust beans) and a chosen fraction of Actigum CS 6 (scleroglucan obtained by aerobic fermentation of a Sclerotium fungus). To approach the diversity of real mucus, seven different gels with scleroglucan concentration from 0.5 to 2wt% were prepared. The mixture was kept under agitation for 48 hours at room temperature. A mass corresponding to 0.2 mL per 10 mL of disodium tetraborate at 0.02 M was added to cross-link the galactomannan chains. The high molecular weight branched macromolecules of sugar, once cross-linked, build a gel matrix that mimicks the mucin network patterning the native mucus. The interested reader can refer to [START_REF] Lafforgue | Thermo-physical properties of synthetic mucus for the study of airway clearance[END_REF], [START_REF] Lafforgue | Rheological properties of synthetic mucus for airway clearance[END_REF], [START_REF] Lafforgue | Rheology of bronchial mucus: characterization and modeling for clearance helping by a medical device[END_REF] for more details.

B. Rheological measurements

The experiments were performed on a stress-controlled rheometer AR 550 (TA Instruments) equipped with a 50 mm -2 • steel cone. Temperature was regulated by a Peltier plate and fixed to 20 • C. Dehydration was prevented by using a wet steel cover insuring a water saturated atmosphere around the sample.

Continuous flow curves as well as three interval thixotropy (3ITT) tests have been performed and the results displayed in [START_REF] Lafforgue | Rheology of bronchial mucus: characterization and modeling for clearance helping by a medical device[END_REF]. To build a full time-dependent model, the flows resulting of multiple structural and load conditions need to be characterized. Hence, series of stepwise changes of shear rate experiments were implemented on samples with 2% of Actigum. Steps at a shear rate chosen as a reference γref are alternated with steps at different shear rates, each step is hold until reaching a steady flow. The method is described in details by Cheng and Evans [START_REF] Cheng | Elastic contributions dominate the viscoelastic properties of sputum from cystic fibrosis patients[END_REF]. These experiments allow to evaluate the thinning and thickening rates on the one hand and the isostructural (constant λ) flow behavior on the other hand. The latter is possible by assuming that the structural response is not instantaneous, this way both the final point of a step and the extrapolation of the subsequent step to the same time provide the response measurement of the same structural state at two different shear levels (see Fig. 1). The experiments were configured so that isostructural flow curves under different reference shear rates were obtained from one single simulant sample, quickly enough to avoid dehydration issues. Flow curves of the sample in its steady structural states at reference shear rates γref from 1.6 (corresponding to patients at rest) to 1000 s -1 (corresponding to cough) were built from the extrapolation of each shear rate step.

III. RESULTS AND DISCUSSION

In the following, all experimental results have been obtained for a 2wt% Actigum concentration, which mimics the best the rheology of cystic fibrosis mucus [START_REF] Lafforgue | Thermo-physical properties of synthetic mucus for the study of airway clearance[END_REF], [START_REF] Lafforgue | Rheological properties of synthetic mucus for airway clearance[END_REF], [START_REF] Lafforgue | Rheology of bronchial mucus: characterization and modeling for clearance helping by a medical device[END_REF]. Fig. 2 displays the continuous flow curve at 20 • C. It confirms the shear-thinning behavior of mucus and the presence of an yield stress (equal to 6.98 Pa). It also highlights a hysteresis loop (±12.5 P a.s 1 ) when performing up and down shear stress ramps.

A. Simplified model

Following the steady-state experiments presented in [START_REF] Lafforgue | Rheological properties of synthetic mucus for airway clearance[END_REF], a Herschel-Bulkley model is chosen as a base to establish a thixotropic model encompassing the time-dependent behavior of mucus:

τ (t) = τ y + K γn (t) (1)
The structure contribution is incorporated by making the rheological parameters depend on a parameter λ(t) that represents the instantaneous structural state. A Herschel-Bulkley behavior is attributed to each structural state:

τ (t) = τ y [λ(t)] + K[λ(t)] γn[λ(t)] (t) + η λ=0 [ γ(t)] γ(t) (2)
τ y is the yield stress, the second term is related to the structural contribution to the viscosity. The term η λ=0 is the residual viscosity when the structure breakdown is at its maximum, this term is currently neglected. The structure kinetic is taken into account from the equation governing the rate of the structure parameter. It consists of three terms describing the competition between breakdown and buildup processes (see equation 3). For the sake of simplicity, a linear dependance on λ and γ is assumed. Such a simplification was found to often provide acceptable approximations according to [START_REF] Mewis | Colloidal Suspension Rheology[END_REF].

dλ dt = -k 1 γλ + k 2 (1 -λ) + k 3 γ(1 -λ) (3)
where k 1 associated to the instantaneous structure λ and shear rate γ traduces the shear-induced breakdown, and k 2 and k 3 associated to the broken structure fraction (1 -λ) (and γ for k 3 ) respectively represent the Brownian buildup and the shearinduced buildup terms.

The kinetic equation for a zero rate of structure evolution leads to the expression of the steady-state structure parameter λ ss :

λ ss = k 2 + k 3 γ (k 1 + k 3 ) γ + k 2 (4)
According to the isostructural curves (not shown here, see in [START_REF] Lafforgue | Rheology of bronchial mucus: characterization and modeling for clearance helping by a medical device[END_REF]), the steady-state isostructural consistency parameter K ss as a function of the reference shear rate shows a rational profile of the form (Fig. 3):

K ss = b + c γ 1 + a γ ( 5
)
with a = 0.0427, b = 13.693 and c = 0.3475. K ss displaying a similar evolution as λ ss , a simple proportional relationship is assumed between these parameters: rate. In other words it quantifies the structuration degree of a stabilized structure under a given shear rate.

K ss = Cλ ss ; λ ss = 1 + k3 k2 γ k1+k3 k2 γ + 1 (6
The value of the apparent yield stress being taken as the extrapolation of the flow curve to zero shear rate, the apparent yield stresses of the isostructure curves discussed above give a trend for the apparent yield stress evolution as a function of the reference shear rate. It is relatively cumbersome to define this relationship due to the lack of data points at high shear rates. From the available data points from 0.029 to 50 s -1 the decrease of yield stress as a function of the reference shear rate is shown in Fig. 5 in addition to the 1000 s -1 point from the alternated SAOS procedure. The other Herschel-Bulkley parameters under steady-state conditions, the Herschel-Bulkley yield stress τ yss and the flow index n ss corresponding to the steady-state structure, can also be expressed as linear functions of λ ss as shown in Fig. 6. Finally the HB parameters are written as follows: The yield stress τ y , the consistency index K and the shearthinning index n are assumed to evolve instantaneously with the level of structure. For instance, the yield stress in a given transitory structural state λ(t) under any shear rate γ1 (t) is equal to the yield stress in the steady structural state λ ss ( γ2 ) corresponding to the same level of structure:

     K ss (λ ss ) =
τ y (λ(t), γ1 ) = τ yss (λ ss ( γ2 )) λ(t) = λ ss ( γ2 ) ⇔ K(λ(t), γ1 ) = K ss (λ ss ( γ2 )) n(λ(t), γ1 ) = n ss (λ ss ( γ2 )) (8)
According to this statement, Eq. ( 2) can be written as follows:

τ (t) = v + uλ(t) + bλ(t) γmλ(t)+p (t) (9)
with v the residual stress at maximum breakdown (λ = 0) (Pa).

The analytical resolution of the kinetic equation ( 3), with t 0 the time at the beginning of the step and λ t=t0 the structural level at this instant, gives:

λ(t) = λ ss + (λ t=t0 -λ ss )e -(k1+k3)γ(t)-k2(t-t0) ( 10
)
t 0 the time at the beginning of the step (s). It can be coupled with the values of C, l 1 and l 2 to only relate λ to γ, t and k 2 :

λ(t) = λ ss + (λ t=t0 -λ ss )e -(aγ(t)+t-t0)k2 (11)
Besides during a given step of the stepwise experiment, γ(t) = γstep = constant, meaning that γ(t) = γstep t and then the expression aγ(t)+t-t 0 within the previous equation becomes:

aγ(t) + t -t 0 = (a γstep + 1) (t -t 0 ) (12)
Thus Eq. ( 11) becomes:

λ(t) = λ ss + (λ t=t0 -λ ss )e -(a γ(t)+1)k2(t-t0) ( 13
)
To get an additional independent equation the expression of λ(t) (13) can be introduced into the equation of τ (t) (9), so that the following system of equations is obtained:

τ (t) = v + uλ(t) + bλ(t) γmλ(t)+p (t) λ(t) = λ ss + (λ t=t0 -λ ss )e -(a γ(t)+1)k2(t-t0) (14)
From Eq. ( 14), the parameter k 2 can be identified to fit the shear stress experimental data. This has been applied to each upwards step from 10 to 50 s -1 with a reference shear rate of 1.6 s -1 . The results are not satisfying because each step is fitted by a quasi linear curve as illustrated by the blue curve of Fig. 7a for the 20 s -1 step. However, the fit becomes very faithful to the experimental data by shifting the τ (λ ss ) curve, changing the v parameter to a value of -0.57 (see the green dash line on Fig. 8 and the red solid line of Fig. 7a). Figure 9 shows that the modification of v (-0.572) offers a good fit of each step (much better than with the initial value v 2 = -3.654). The 10 -1 step is not fitted as well as the others with a maximum error around 5%.

In order to obtain a better fit of each curve, a non-constant value of v (denoted as v 2 ) would need to be taken into account as a function of the reference shear rate. Tab. I displays the variations of v 2 versus the reference shear rate.

At this stage, it seems thus that the time dependent rheological state equation assumed in Eq. ( 9) is not sufficient to account for the whole set of experimental data.

B. Full model

A more complete model is proposed (see Eq. ( 15)) in the present section. It accounts for the viscoplasticity at Fig. 9. Shear stress experimental data (circles) and Eq. ( 14) with v changed to -0.572 (solid lines). extreme structural degrees, respectively τ y λ=1 corresponding to the yield stress of the native structure (λ = 1) and τ y λ=0 corresponding to the yield stress of the entirely broken down structure (λ = 0). The shear-thinning behavior is modeled by the power-law term and the residual viscosity η λ=0 when the breakdown of the material structure is maximum (λ = 0) is accounted for through the last term of the first equation.

In comparison to Eq. ( 11), the expression of the time evolution of the structure degree λ suffers here of no more simplification and is governed by the three kinetic constants k 1 , k 2 and k 3 . The complete model writes: k1+k3) γ(t)+k2](t-t0) (15) A Matlab program involving a global optimization algorithm was implemented to fit the latter model to several sets of experimental data. The experimental data to fit are stress responses to series of stepwise changes of shear rate from 1.6 s -1 (in the range of resting respiratory conditions) to 100 s -1 (in the range of coughing conditions). These steps were implemented from reference shear rates of 1.6, 10, 20, and 50 s -1 . This means that upward and backward steps (respectively destructuration and restructuration steps) are involved here from various structural states of reference. The algorithm was designed to minimize the objective function F obj (Eq. ( 16)) involving the square of the differences between each p to N p data point τr,s,p , of each s to N s step, of each r to N r set (corresponding to each reference shear rate) and the corresponding model prediction τ r,s,p : The experimental data and the corresponding model predictions are presented in Fig. 10. This result was obtained by minimizing the error F obj to the best value of 0.936 corresponding to the parameter values gathered in Eq. (17).

     τ (t) = (τ y λ=1 -τ y λ=0 ) λ(t) + τ y λ=0 + Cλ(t) γmλ(t)+p (t) + η λ=0 γ(t) λ(t) = λ ss + (λ t=t0 -λ ss )e -[(
              τ y [λ(t)] = (τ y λ=1 -τ y λ=0 ) λ(t) + τ y λ=0 K [λ(t)] = Cλ(t) n [λ(t)] = mλ(t) + p η λ=0 = 0.003 λ(t) = λ ss + (λ t=t0 -λ ss )e -[(k1+k3) γ(t)+k2](t-t0) (17)
with τ y λ=1 = 20.45, τ y λ=0 = 5.30, C = 6.345, m = 0, p = 0.345, k 1 = 0.0012, k 2 = 0.2542 and k 3 = 0.0014.

The modified Herschel-Bulkley model proposed here is thus robust enough to predict faithfully numerous sets of buildup and break-down transient data in terms of viscoplasticity, shear-thinning and thixotropy. The model prediction of the corresponding time evolution of the structure degree λ is presented in Fig. 11. It shows values contained between 0 and 1 with increasing values toward a plateau for downward shear rate steps (build-up) and decreasing values toward a plateau for upward shear rate steps (break-down).

IV. CONCLUSION

New rheological measurements using an hybrid rheometer have been proposed for a synthetic mucus simulant composed of Actigum and Viscogum. At a concentration of 2wt% in Actigum, the polymer mimics human sputum from cystic fibrosis patient. Mucus is a highly complex material exhibiting, from simple flow curve, an yield stress, shear thinning, viscoplastic and thixotropic properties. Series of stepwise changes of shear rate experiments have been also performed to obtain the steady-state flow curves corresponding to several isostructural states that were fitted with Herschel-Bulkley models. A quantification of the structure parameter corresponding to the stable structural states was done and the structure dependency of the rheological parameters was identified. The model has been then improved to account for the thixotropic structure kinetic (time-dependent rheological state equation). The simplified model showed in average, a pretty good fit of the transitory steps up from a reference at rest taken as 1.6 s -1 (respiration conditions) up to 100 s -1 (cough) shear rate conditions. The proposed improved timedependent model accounts for the existence of a shear-rate dependent supplementary stress related to a residual viscosity for maximum break-down. This model is robust enough to predict faithfully numerous sets of build-up and break-down transient data in terms of viscoplasticity, shear-thinning and thixotropy.

An extension of this work is to further validate the full model by new experiments on real bronchial mucus collected from patient. Another perspective is to use the present model in numerical simulations of the transport of bronchial mucus in an idealized bronchus, to account for the time-dependent nature of the mucus.

Fig. 1 .

 1 Fig. 1. Profile of a stepwise change in shear rate.

Fig. 2 .

 2 Fig. 2. Typical rheogram for a 2wt% concentration in Actigum obtained at 20 • C. Continuous shear stress ramp (up and down).

Fig. 3 .

 3 Fig. 3. Evolution of the consistency index Kss as a function of the reference shear rate.

  ) with C a constant. One can then identify and assign the values for: C = b = 13.693, l 1 = k 3 /k 2 = c/b = 0.025 and l 2 = k 1 /k 2 = a -c/b = 0.017. Fig. 4 shows the evolution of the steady-state structural parameter as a function of the shear

Fig. 4 .

 4 Fig. 4. Evolution of the steady state structure parameter as a function of the reference shear rate for a 2% Actigum sample.

Fig. 5 .

 5 Fig. 5. Evolution of the apparent yield stress as a function of the reference shear rate for a 2% Actigum sample. The values from 0.029 to 50 are drawn from the stepwise experiments.

Fig. 6 .

 6 Fig. 6. Evolutions of (a) τyss and (b) nss as functions of λss.

Fig. 7 .

 7 Fig. 7. (a) Fit of Eq.(14) with the shear stress experimental data (blue solid line) and fit of the same equation with v changed to -0.572 (red solid line); (b) Corresponding evolutions of the structural parameter.

Fig. 8 .

 8 Fig. 8. Impact of the modification of v on the τ (λss) curve (green dash line).

  ,s,p -τ r,s,p )

Fig. 10 .

 10 Fig. 10. Shear stress experimental data (symbols) in response to stepwise changes of shear rate and corresponding model predictions from Eq. (15) (solid lines). Each graph corresponds to a reference shear rate from 1.6 to 50 s -1 respectively from the left to the right.



Fig. 11 .

 11 Fig. 11. Model predictions of the evolution of the structure degree as a function of time from Eq. (15). Each graph corresponds to a reference shear rate from 1.6 to 50 s -1 respectively from the left to the right.

TABLE I IDEAL

 I VALUE OF v FOR EACH STEP: v 2 ( γref )

	γref (s -1 )	v 2 (Pa)
	10	-1.6296
	20	-0.5778
	30	-0.2962
	40	-0.1334
	50	-0.2231

ACKNOWLEDGMENT

The authors are indebted to the Laserson Co. for kindly providing the Actigum and Viscogum polymers.