Using rheological monitoring to determine the gelation kinetics of chitosan-based systems

1.   Introduction

The characterization and application of polymers in Biotechnology, Biomedicine and Therapeutics has expanded exponentially over the last decades. Particularly, biopolymers have been widely and increasingly studied and employed in broad range of applications. They became of great interest in some growing areas of knowledge such as tissue engineering and nanotechnology, since the intrinsic structure of polymeric networks provides interesting properties that can be practically exploited ^[2,3,4,5,6]. Chitosan is a typical biopolymer that falls in this class.

Chemically, chitosan is a primary aliphatic amine obtained by alkaline deacetylation of chitin found in skeletal structure of crustaceans, insects, mushrooms as well as in cell wall of fungi ^[7]. The characteristics of chitosan are deeply analyzed in Section 3 but the main properties influencing its characteristics (see the chemical representation in Figure 1(a)) are the average molecular weight (that measures the average size of the chain) and the degree of deacetylation (that represents the proportion of deacetilated units in matrix) ^[8]. The pH substantially alters the charged state of the monomer units, and consequently the properties of the chain. At low pH (below 6) the amino groups are positively charged, and the chitosan is a water-soluble cationic polyectrolite ^[8]. Contrarily, at higher pH (above 6.5) the amino groups become deprotonated; consequently, the polymer loses its charge and becomes water-insoluble. When solubilized the repulsion forces are weak, which allows the inter-association between polymer units, that can yield fibers, films or hydrogels, depending on the conditions employed to initiate the modulation. The reactivity of amino groups enables a broad range of reactions that can be employed for functionalization and/or crosslinking with the purpose of conferring elasticity and chemical stability, which enhance its use in tissue engineering applications and drug release matrices ^[9,10,11]. Nevertheless, the application of chitosan based matrices has been limited due to its poor biomechanical properties and uncontrolled degradation rate in vivo. Several modifications of the nature of the linkages between the polymeric units have been attempted with the purpose to improve the structural properties of the matrices, especially when they are in gel form.

It is now accepted that chitosan is capable of forming thermally stable gels in presence of crosslinking agents or other negatively charged molecules. Several authors investigated the use of genipin for the purpose of controlling gel properties through the manipulation of crosslinking mechanisms ^[12,13,14]. Genipin is a naturally occurring crosslinking agent that induces the formation of crosslinking strands between chitosan molecules, thus improving the mechanical properties of the gel formed. Furthermore, it is capable of binding with other biopolymers such as gelatin, and able to be used enclosed in drug delivery matrices since it has a lower toxicity and good degradation properties ^[14].

In order to achieve a better understanding of the gelation mechanisms and derive models that can be useful for forecasting and control purposes, several works have been reported aiming to characterize pure chitosan self-crosslinking kinetics and agent-based crosslinking (gelation) kinetics, see Merkovich et al. (2001) ^[15] and Vo et al. (2020) ^[16] among others. Rheological monitoring is an analytical technique commonly used to infer about the polymer reaction kinetics, since it can be carried out in situ ^{[17,18,19,20,21]}. The usefulness of rheological behavior to study polymerization processes and construct kinetic models is enormous, see Thévenot et al. (2007) ^[22] among others. However, its potential is not fully exploited because fitting theoretically or empirically based models to rheological data for parametrization purposes is still performed based on linear regression procedures combined with logarithmic transformations ^[23,24].

The parameter estimation of models comprising sets of Differential-Algebraic Equations (DAEs) is now state of the art for process systems engineering community. A vast range of applications can be found in the open literature (see Tjoa and Biegler (1991) ^[25] and Samaniego et al. (2007) ^[26] among others). The application to fit reaction kinetics employing the integral method ^[27,28] combined with robust numerical procedures standing on optimization ^[29] is extensively found in chemical reaction literature but still scarce when applied to polymeric systems. Aiming to fulfill the gap regarding the parametrization of kinetic rates of crosslinking reactions we apply the approach to models representing the dynamics of gelation of chitosan systems using rheological data.

Several authors found evidences that the rheological behavior of the polymers involved in crosslinking reactions is highly correlated with the rheological degree of conversion of liquid phase to gel phase (or solid phase in other cases). This relation was exploited to interpret the kinetics of crosslinking reactions of polymers ^[30], and to relate the rheological degree of conversion with the reaction kinetics. Here, we also exploit it for monitoring the kinetics of self-crosslinking and agent-based crosslinking of chitosan. Several works parametrize the kinetics of hybrid crosslinking of chitosan with genipin using the gel time, that corresponds to the time instant at which the elastic and viscous moduli become equal. Then, considering a set of at least two distinct temperatures the activation kinetic constant and the reaction order can be determined from the equations proposed by Butler et al. (2003) ^[31] or Ross-Murphy (1991) ^[32]. This approach provides useful estimates but has three main drawbacks: ⅰ. it requires characterizing the rheological behavior of the system at temperatures different of 37℃ which is unrealistic for some applications; ⅱ. the time at which the gel point is attained is sometimes hard to detect; and ⅲ. the parameter estimates are based on a limited number of data points corresponding to the temperatures tested. We propose and demonstrated a strategy that overcomes the main issues listed as it infers the kinetics from a complete dynamic rheological test (at a pre-specified temperature) which increases the accuracy of the parameters obtained, and avoids the error contamination of parameter estimates based on tests at distinct temperatures.

The paper includes three innovative elements: i. the application of systematic approaches to fit kinetic models of gelation reactions using the integral method with chitosan gels being used as demonstrating examples; ii. the use of rheological data to fit (empirical) kinetics in terms of the rheological degree of conversion; and iii. the use of models of the family of ^[1] postulated from rheokinetics studies for describing the transition dynamically.

In the remaining sections, bold face lowercase letters represent vectors, bold face capital letters are for continuous domains, blackboard bold capital letters are for discrete domains and capital letters are for matrices. Finite sets containing $\iota$ elements are compactly represented by $\left[\!\left[ \right.\!\right. \iota \left.\!\left. \right]\!\right] \equiv\{1, \cdots, \iota\}$. The transpose operation of a matrix/vector is represented by "$^\intercal$", and the circumflex is used to designate expectation.

This paper is divided into four additional sections. In Section 2 the parameter estimation algorithm employed is presented. Section 3 introduces the materials and the preparation procedure used for producing and characterizing gels, and overviews the fundamentals behind the reokinetics. Section 4 presents the results of model fitting of the reaction rate for self-crosslinking gelation of chitosan and agent-based crosslinking of chitosan using genipin. Finally, Section 5 presents the conclusions.

2.   Parameter estimation in dynamic models

This section overviews the parameter estimation of dynamic models (i.e., represented by evolutive differential equations). First, we introduce the model fitting problem using Maximum Likelihood criterion. Then, we establish the basis to numerically construct confidence intervals for parameters, analyze the parametric collinearity and check model identifiability.

The parameter estimation algorithm employed considers general processes described by DAEs systems with the form:

where ${\bf{f}}(\bullet)$ and $\bf{g}(\bullet)$ are vectors of continuous functions, ${\bf{x}}(t)\in\mathbb{R}^{n_s}$ is the vector of state variables (which characterize the dynamics of the system), $\dot{{\bf{x}}}(t)$ their time derivatives, ${\bf{y}}(t)\in \mathbb{R}^{n_o}$ is the vector of variables measured (or observed), ${\bf{u}}(t) \in \mathbb{R}^{n_i}$ the vector of input variables, ${\bf{p}} \in \mathbb{R}^{n_p}$ is the set of parameters in the phenomenological model (2.1a), $\pmb{\theta} \in \mathbb{R}^{n_\theta}$ is the complete set of parameters to fit which includes ${\bf{p}}$ and those in the measurement error variance model, ${\bf{x}}_0$ the vector representing the initial state of the system and $\pmb{\epsilon}$ the measurement error. Further, $n_s$, $n_o$, $n_i$, $n_p$, and $n_\theta$ are the number of states, outputs, inputs, parameters in the model, and of the complete set of parameters, respectively. The model represented by the set of Eq (2.1) is a nonlinear state space model where Eq (2.1a) is the state equation, Eq (2.1b) the output equation and Eq (2.1c) the initial condition.

Here, we consider that the objective function of the model parametrization problem is the Maximum Likelihood Estimate (MLE) criterion. Several applications found in literature use the Least Squares criterion (see for example Faber et al. (2003) ^[33]), however the MLE is preferable for non-normal data or when the probability distribution function (pdf) of observations is unknown or difficult to compute ^[34]. Rheological data are often affected by varying observational error and the MLE is expected to assure better parameter estimates and simultaneously characterize the error model underlying the measurement system. When the observations are independent on each other and normally distributed with constant variance both criteria have the same performance. For such a situation the maximization of log-likelihood is equivalent to minimize the residual sum of squares. In this work we assume that the measurement error is normally independent. However, to cope with the differences of range of the rheological properties measured during the experiments, where final values easily achieve three orders of magnitude above the values at the beginning, an heteroscedastic model of the variance of the forecasts is considered. That is, the increasing trend of the prediction error is accounted for in model fitting and we consider that the variance of each variable, $\sigma^2$, follows the general heteroscedastic model

where $w$ and $\gamma$ are additional parameters to fit. Let $\pmb{\beta}\equiv\{w, \gamma\}$; then $\pmb{\theta}\equiv{\bf{p}}\cup\pmb{\beta}$, $\mathbb{R}^{n_\theta}\equiv \mathbb{R}^{n_p+2}$ and $n_\theta = n_p+2$.

For the sake of generalization let us consider a set of variables measured along the experiments, containing various outputs and inputs, represented by vector ${\bf{z}}(t)\in \mathbb{R}^{n_v}$. Specifically, ${\bf{z}}(t)\equiv\{{\bf{y}}^m(t)\in \mathbb{R}^{n_{o, m}}, {\bf{u}}^m(t)\in \mathbb{R}^{n_{i, m}}\}$, with $\mathbb{R}^{n_{o, m}} \subseteq \mathbb{R}^{n_o}$ being the domain of observed variables, $\mathbb{R}^{n_{i, m}} \subseteq \mathbb{R}^{n_i}$ the domain of input variables observed; the subscript "$m$" stands for measurement. Further, $n_{o, m}$ is the number of output variables measured, $n_{i, m}$ the number of inputs measured in the experiment over the time, and $n_v = n_{o, m}+n_{i, m}$ the number of variables measured. The measurements are designated as $z_{i, j}$ with the subscript $i$ standing for the variable measured after convenient ordering, and $j$ for the discrete time instant at which the process is sampled. Practically, the measurements are obtained by sampling/observing the system at discrete instants $t_j, j\in \{0, \cdots, n_t\}$, with $n_t+1$ being the number of time instants at which data is available. The measurement error of each variable, designated as $e_{i, j}$, is assumed to be normally independent with variance $\sigma^2_{i, j}$. Consequently, the log-likelihood function is

To distinguish between values predicted by the model and measurements used for model fitting, the former are represented by vector $\hat{{\bf{z}}}$, and the latest by ${\bf{z}}$. Therefore, the optimization problem to solve is as follows,

where Equation (2.4a) is the log-likelihood function, Eq (2.4b) represents the state space model to fit, Eq (2.4c) constructs the set of model predictions used for fitting, Eq (2.4d) is the set of measurements available, Eq (2.4e) is the variance model, Eq (2.4f) bounds the parameters, $N$ is the number of measurements available, $\pmb{\theta}^L$ the parameters lower bound, $\pmb{\theta}^U$ the upper bound, and $\varepsilon$ a small positive constant chosen to avoid numerical problems. When $\gamma = 0$ in Eq (2.4e), the model falls into a constant variance form; contrarily, when $\gamma = 1$ it reduces to constant relative variance.

There are two distinct strategies for addressing the parameter estimation in dynamic models. The first uses integration codes to determine the dynamics of observed variables and the sensitivities with respect to the set of parameters. The model formed by Eq (2.1) is solved simultaneously with four additional equations that enable computing the parametric sensitivities (derivatives of observed variables with respect to parameters, i.e., $\rm{d}{\bf{y}}/\rm{d}\pmb{\theta}$). We note that the sensitivity is represented by a $n_v\times n_\theta$ matrix, and the equations to solve are as follows

where $B({\bf{x}}, \pmb{\theta}, t)$ is the matrix of sensitivities of the state variables with respect to ${\bf{p}}$, $S({\bf{x}}, \pmb{\theta}, t)$ the sensitivity of the observed variables with respect to the overall set of parameters, $\mathcal{M}({\bf{x}}, \pmb{\theta}, t)$ the Fisher Information Matrix (FIM) and $C$ the measurement error covariance matrix. Model predictions and parametric sensitivities are employed to build the objective function, gradient and Hessian matrices required by optimization tools. Non-Linear Programming (NLP) solvers are required to handle the problem formed by Eqs (2.4) and (2.5), see Stewart et al. (1992) ^[35]. This approach is generally designated as sequential because it sequentially iterates between DAEs model solution and likelihood minimization until convergence.

The second group of methods, called simultaneous-based approaches, employ discretization techniques, such as orthogonal collocation on finite elements, to transform the original model into a set of algebraic equations, one per experiment and measurement available. The problem is then solved via a NLP code in the sense of Maximum Likelihood or Least Squares ^[36,37]. This approach leads to a dynamic optimization control problem where the goal is to choose the parameter set that minimizes the residual sum of squares.

Here, we employ the sequential approach since it is easier to code and successful for large scale NLP problems ^[38,39]. The DAEs system is solved via a Backward Differentiation Formula (BDF) based algorithm - DASPK3 ^[40]. This solver has the ability of determining simultaneously the model predictions and the sensitivities using Eqs (2.1), (2.2) and (2.5). We note that the jacobian of the optimization problem can be found by analytical differentiation of Eq (2.1). The optimization code employed to handle the optimization problem was the bound-constrained quasi-Newton method L-BFGS-B ^[41].

After solving the optimization problem formed by Eqs (2.4) and (2.5) various metrics of quality of fit are determined. Together, they intend to assess the model accuracy and model adequacy to represent the underlying data. Typically, this post-analysis is performed based on FIM, which has important limitations, especially for nonlinear models where it only provides a lower bound for variance and symmetric parametric confidence intervals ^[34]. However, alternative approaches such as bootstraping ^[42], despite of the advantages, are computationally intensive, especially when the amount of data is large.

To compute $100\times (1-\alpha)\; \%$ confidence level intervals for parameters we use F probability distribution function to provide upper bounds for square differences between real values ($\pmb{\theta}$) and estimates ($\hat{\pmb{\theta}}$), i.e.,

where $F(\alpha, n_{\theta}, n_t+1-n_{\theta})$ is the value of the $F$ distribution with $(n_{\theta}, n_t+1-n_{\theta})$ degrees of freedom at $1-\alpha$ and $s$ is an approximation of the standard error. Practically, it is given by

The computation of the variance-covariance and correlation matrices of the parameters is based on the procedures in ^[43,44]. In practice, we construct the parameter variance-covariance matrix by inverting the FIM; that is

Next, we use $\mathcal{C}({\bf{y}}, \pmb{\theta}, t)$ to determine the parametric correlation matrix,

where $\operatorname{diag}(\mathcal{C}({\bf{y}}, \pmb{\theta}, t))$ is the matrix formed by the diagonal elements of $\mathcal{C}({\bf{y}}, \pmb{\theta}, t)$. The local identifiability of the model is checked using information of i. the FIM ^[45]; and ii. parameter variance-covariance matrix. We first determine the minimum eigenvalue of the FIM, $\lambda_{\min}[\mathcal{M}({\bf{y}}, \pmb{\theta}, t)]$, and in case it is larger than 1 × 10⁻⁵ the system is considered practically identifiable. The occurrence of $|\mathcal{R}_{i, j}|\geq 0.9$, where $\mathcal{R}_{i, j}\in\mathcal{R}({\bf{y}}, \pmb{\theta}, t), \; i\neq j$ is the correlation between parameters $i$ and $j$, is an indication of collinearity between both parameters ^[46].

Finally, we check the model adequacy by using the $\chi^2$ goodness of fit test. The absolute and relative tolerances imposed to integrator are 1 × 10⁻⁸ and 1 × 10⁻⁷, respectively. The absolute and relative tolerances imposed to optimizer are 1 × 10⁻⁶ and 1 × 10⁻⁵, respectively.

3.   Kinetics of the gelation reaction of pure chitosan solutions

This section introduces the theory supporting the formation of hydrogels and the background behind the use of rheological information for estimating the reaction rate in polymer-based transformations. In Section 3.1 we address the formation of hydrogels and Section 3.2 overviews the materials used in experiments and sample preparation procedure. Further, we also characterize the device used for rheological monitoring.

We note that chitosan is a copolymer containing $\beta$-1, 4-N-acetylglucosamine, see the schematic representation in Figure 1(a) ^[7].

3.1.   Gel formation mechanisms and rheological characterization

Here, we analyze the mechanisms supporting the transformation of (bio)-polymer-based materials from a dispersion to a gel phase via crosslinking, and the use of rheometry for monitoring the phenomena.

Hydrogels based on cross-linked chitosan can be divided into three types regarding their formation mechanisms: i. self-crosslinking; ii. hybrid polymer networks commonly originated by using a cross-linker; and iii. interpenetrating polymer networks ^[47]. The first results from the combination of two structural units of chitosan which may (or may not) belong to the same network. In hydrogels based on hybrid polymer network the crosslinking reaction occurs between a structural unit of the chitosan chain and a structural unit of a polymeric chain of another polymer. Finally, interpenetrating polymer networks contain a non-reacting polymer added to the chitosan solution before crosslinking, thus producing a cross-linked network in which the non-reacting polymer is entrapped. In Section 4.1, we investigate the self-crosslinking reaction, and in Section 4.2 the study addresses the kinetics of crosslinking in hybrid polymer network formation employing the cross-linker agent genipin.

The gel formed by chitosan reaction with a crosslinking agent, here represented by genipin, which structure is in Figure 1(b), is widely studied due to its application in biomedical and food engineering areas. The gelation mechanism was identified by ^[48], and involves the formation of hybrid cross-linked networks of heterocyclic entities with 1 to 4 amine units at physiological pH.

The study of the kinetics of ordinary chemical transformations involving low molecular weight reactants and products generally lead to small variations of the viscosity of the reaction mixture. The viscosity depends not on the degree of conversion, only on the reaction conditions at which it is carried out. Contrarily, polymerization and depolymerization reactions involve the formation and disappearance of high molecular weight components. Thus, their kinetics will have a strong effect on the viscosity of the reaction mixture. The study of reaction kinetics based on the monitoring of the rheological properties of the reaction mixture gave origin to a new field of knowledge designated as rheokinetics ^[1].

The reaction rate is generally studied by measuring the concentrations of reactant species employing techniques such as spectrophotometry, differential calorimetry, and chromatography. Rheokinetics uses the changes in viscoelastic properties to estimate kinetic parameters and to validate the reaction mechanisms for mixtures of polymers. Oscillatory deformation tests are easily applied in rheometry equipment allowing to follow the mixture behavior and study the underlying kinetics ^[49]. The shear modulus, represented as $G^{*}$, is resolved into its real, $G^{'}$, and imaginary, $G^{''}$, components:

where $G^{'}$ is the elastic modulus, which measures the energy stored per oscillation cycle, and $G^{''}$ is the viscous modulus that represents a measure of the energy dissipated as heat per oscillation cycle. Hydrogels formed by crosslinking reactions denote a rheological behavior characterized by the dominance of the elastic modulus for extended reaction time. This behavior is frequency independent at the gel point, identified as the instant at which the moduli are equal ^[50].

Rheokinetics is commonly used for monitoring the kinetics of polymerization reaction; first it was applied to monitoring resins curing process ^[51] and later extended to other systems, including systems involving sol-gel transition similar to that addressed in this paper, see ^[52]. The rheokinetics allows monitoring of viscosity dynamics and the elastic modulus reflects the characteristics of this behavior. Thus, we omit the information relative to viscous module as it is irrelevant for kinetics parametrization, typically denoting a constant behavior after the gel point.

3.2.   Material, sample preparation and rheological measurement

Here we present the experimental procedure employed in monitoring the rheological behavior of chitosan solutions that produce data used for determining the kinetic rate of transformation of the reactants into gel. The chitosan used in the experiments was purchased from Sigma-Aldrich in powder form, with a molecular weight of approximately 2 × 10⁵ Dalton and a degree of deacetylation greater than 85 %. Hydrated $\beta$-glycerol-phosphate disodium salt (C₃H₇Na₂O₆P · xH₂O; FW = 218.05 g/mol g) used to adjust the pH of the chitosan solutions was also purchased from Sigma-Aldrich. Genipin (crystal-like powders, reagent grade) was supplied by Challenge Bioproducts Co., Taiwan. The chitosan hydrogels considered in this study are to be prospectively used in tissue engineering and as drug release matrices. The $\beta$-glycerol-phosphate disodium salt is commonly used to modulate the pH of chitosan-based systems gelation to the human body range, i.e., 7.35–7.45, see ^{[53,54,55,56,57]} among others. In our study, this salt was also used for the same purpose. All other reagents and solvents used were of analytical grade.

First, an aqueous solution of chitosan was prepared by dissolving 2 g of chitosan in 100 cm³ of distilled water containing acetic acid (0.5 % (v/v)), at room temperature. The chitosan solution was then filtered and subsequently sonicated. Thereafter, a given amount of glycerol-phosphate disodium salt was dissolved in distilled water and carefully added drop by drop, under magnetic stirring, to the chitosan solution to obtain a clear and homogeneous liquid with concentration of 1.5 g chitosan/100 cm³ and pH of 7.0. Genipin powder was dissolved in the resulting solution to produce samples of genipin concentration of 0.15 % (w/w). Chitosan solutions without genipin were used as control.

The rheological characterization of the hydrogel samples was performed in a C-VOR rheometer from Bohlin Instruments, Inc., USA equipped with a cone-and-plate geometry (cone angle of 4° and diameter of 20 or 40 mm) provided with a circulating environmental system for temperature control. The rheometer used is a shear stress controlled instrument. That is, it applies a controlled shear stress (torque) and measures the resultant displacement represented as the rotation rate. Torque and movement are converted into "rheological format" (i.e., $G^{'}$ vs. $G^{''}$) using the calibration system, and exported to MS Excel. To prevent drying of the samples during experiments, a steel ring of approximately 60 mm is placed around the measuring geometry, the annulus was filled with water, and the sample-holding region was sealed with an insulated cover. As soon as the sample was introduced onto the plate data started being gathered with this time instant being the reference for subsequent treatment. A time sweep in the linear viscoelastic region, at a constant shear frequency of $\omega$ = 1 Hz, was performed to monitoring the in situ gelation behavior of the hydrogel dispersions at constant temperature of 37 ℃ and physiological pH. The stress load employed is 1 Pa with the elastic and viscous moduli dynamics being monitored over the time of the experiment.

4.   Model fitting of gelation kinetics

In this section we apply the methodology described in Section 2 to fit the kinetics of chitosan crosslinking using rheological data. In Section 4.1, we consider the self-crosslinking mechanism, and in Section 4.2 we consider the hybrid network formation mechanism resulting from using genipin as cross-linker.

The model used to describe the chitosan gelation reaction monitored rheologically is owned to ^[1], and has been often applied to describe the dynamic behavior of similar polymerization systems (see for instance Zlatanic et al. (1999) ^[58] and Madbouly and Ougizawa (2004) ^[59]). Its structure is purely empirical, and is employed to explain the rheological behavior of polymer reactions occurring by crosslinking of monomer units of the same network or hybrid polymerization using units of different networks where an external agent plays the role of cross-linker. The model takes the form

where $\alpha$ is the rheological degree of conversion of chitosan into gel, $n$, $k_1$ and $k_2$ are parameters, with the term $n+1$ providing an indication of the reaction order. Eq (4.1a) represents the kinetic rate, and Eq (4.1b) is for the initial condition. The model represented by Eq (4.1) is employed to establish relations between the rheological behavior of the polymers and the chemical transformations occurring, even in early stages of the reaction. In radical polymerization reactions $\alpha$ is proportional to the degree of conversion of the monomer. Similarly, for gelation reactions $\alpha$ stands for the extent of polymer that turns into gel phase ^[60]. The rheological degree of conversion is to be evaluated from the elastic modulus employing the relation

where $G_{\rm{scl}}^{'}(t)$ is the elastic (storage) modulus at instant $t$, $G_{\rm{scl}}^{'}(0)$ its value at the beginning of the experiment, $G_{\rm{scl}}^{'}(t_f)$ its value at the time the experiment ends ($t_f$), and the subscript "$\rm{scl}$" denotes self-crosslinking.

The experiments employed to fit the model formed by Eqs (4.1) and (4.2) lasted for about 650 min (i.e., $t_f$ = 650 min) comprising 1321 observations (the frequency of data acquisition is 30 sec). Several works refer that the gelation process of polymers evolves into three sequential phases: ⅰ. the initiation, where the elastic modulus remains close to 0; ⅱ. the sol-gel transition, characterized by large increases of the elastic modulus; and ⅲ. the equilibrium, where the network almost achieved its steady-state and there is a slight increase of $G^{'}$ ^[17,61]. A similar behavior was observed for the gelation of chitosan and chitosan cross-linked with genipin dispersions. Particularly, the first two phases are well identified, see Figure 2 where the data is plotted, but the third can not be identified for this experimental horizon.

We note that the gelation curves present a discontinuity of the value of $G^{'}$ followed by a nearly steady state where the equilibrium is achieved. The behavior is observed at around 650 min for self-crosslinking of chitosan and at 135 min for hybrid cross-linked chitosan solution. Various authors provide evidences that it can be due to the degradation of the sample structure by synerysis, a phenomena already identified by ^[62] for chitosan-based hydrogels. Since the rheological behavior of the solutions beyond the discontinuity is mainly due to the mechanical phenomena involved, only the measurements gathered before were used to fit the kinetic rate. This strategy discards the phase during which the shrinkage of the gel and water expelling do occur with induced mechanical transformations dominating the chemical phenomena ^[63].

4.1.   Model of the self-crosslinking gelation reaction of chitosan

Here, we fit the kinetic rate for self-crosslinking of chitosan. Table 1 presents the values of the parameters and corresponding 95 % confidence intervals obtained with the algorithm described in Section 2. The reaction order found for chitosan self-crosslinking kinetics is $n+1 = 1.3375$. This value is smaller than that determined by ^[59] for poly(vinyl methyl ether) crosslinking employing rheological data and a kinetic model with similar structure to ours, and by ^[64] for n-isobutyryl chitosan hydrogels employing a model based on the gelation time ^[32]. The order of the kinetic rates found in both works is 2. Figure 3 compares the model predictions with experimental data and demonstrates the model accuracy in explaining the behavior of the self-crosslinking reaction of chitosan in the first two phases of the gelation process.

The model initially proposed by ^[1] to describe polymer reactions via rheokinetics (formed by Eqs (4.1) and (4.2)) includes an additional parameter - the power of term $(k_1+k_2 \ \alpha)$. In this study, we fixed it to 1 to avoid parameter collinearity, and because from phenomenological point of view it is straightforward to consider that the increase of the extension of the dispersion that turn into gel state over the time is proportional to the amount that already is in gel phase, see Mi et al. (2005) ^[48]. To analyze our assumption we used the $\chi^2$ goodness of fit test, which proved the model adequacy in describing the underlying process. In practice, considering an additional parameter might lead to local identifiability issues. The three parameter-based model here fitted is fully observable ^[34] and locally identifiable using the criterion introduced in Section 2. We obtained $\lambda_{\min}[\mathcal{M}({\bf{y}}, \pmb{\theta}, t)] =$ 2.545 × 10⁻² well above of the threshold value imposed.

The lower triangular part of parameters correlation matrix, $\mathcal{R}({\bf{y}}, \pmb{\theta}, t)$, is shown in Table 2, and reveals the linear independence of the parameters on each other, thus strengthen the conclusion about the local identifiability of the model. If an additional parameter is included in the model, such as the above referred power index, two or more of the parameters would become auto-correlated; consequently, the quality of the fitting does not improve significantly but the model parsimony degrades ^[65] as well as the model performance in forecasting scenarios.

The parameter estimation procedure also reveals that the optimal values of $w$ and $\gamma$ are 0.04089 and 1, respectively, providing evidence that the variance of the prediction is optimally represented by a constant relative variance model. Finally, this result demonstrates the advantages of using the MLE criterion for fitting rheological data as it allows obtaining the best model for the prediction error variance accounting for eventual trends in measurement error pdf.

4.2.   Kinetics of hybrid network formation

Now, we fit the kinetics of the gelation of chitosan-genipin solution 0.15 % (w/w). We are interested in discriminating the self-crosslinking kinetics from the hybrid crosslinking induced by chitosan. Thus, we adapt the Eqs (4.1) and (4.2) considering that two terms contribute to the behavior of the elastic modulus. The first is due to self-crosslinking, and is represented by the difference $(G_{\rm{scl}}^{'}(t)-G_{\rm{scl}}^{'}(0))$. This component is determined from pure chitosan experiments and was considered in Section 4.1 to fit self-crosslinking kinetic rate. The second term is the contribution of genipin via hybrid crosslinking represented by the difference $(G_{\rm{hcl}}^{'}(t)-G_{\rm{scl}}^{'}(t))$. We note that subscript "hcl" designates hybrid crosslinking and without loss of generality we assume that $G_{\rm{hcl}}^{'}(0) = G_{\rm{scl}}^{'}(0)$, as shown in Figure 3. Herein, the rheological degree of conversion, represented by $\beta$, is

with $t_f =$135 min (see Figure 3). The equation representing the rheological degree of conversion due to hybrid crosslinking is

The model including Eqs (4.3) and (4.4) was fitted with the algorithm described in Section 2. Table 3 presents the values of the parameters fitted with $w$ and $\gamma$ again indicating that the data is also optimally fitted by constant relative variance models. The order of the hybrid crosslinking reaction is $n+1 = 2.2402$, a value that is near of the average of chitosan units linked per heterocyclic unit ^[48]. The reaction order is larger than the value of 1.78 obtained by ^[31] employing the reciprocal of gel time coupled with the equation proposed by ^[66]. However, the former authors refer that the value obtained is lower than that expected from an irreversible gelation reaction (order 2) due to the occurrence of secondary reactions. Our value of reaction order almost doubles the value achieved for pure chitosan gelation kinetics, thus highlighting the advantages of heterogeneous crosslinking in reducing the gelation time, as was demonstrated by ^[18]. Dimida et al. (2015) ^[21] obtained $n\in[1.83, \; 1.90]$ by fitting gel time information obtained at different temperatures. Using a similar approach Espinosa-García et al. (2007) ^[19] obtained $n = 3.02$. Finally, Fraga et al. (2007) ^[67] uses Differential Scanning Calorimetry data combined with a Least Squares-based fitting approach, and obtained $n = 2$. This last study considers the model of Sourour and Kamal (1975) ^[68] for fitting which has structural similarities with that of Malkin and Kulichikin (1996) ^[1] used in our work. This resemblance is particularly interesting as it allows comparing values obtained from structurally comparable models. All these values are in good agreement with those obtained herein.

The model adequacy was also confirmed with the $\chi^2$ goodness of fit test. To infer about the model identifiability we determined $\lambda_{\min}[\mathcal{M}({\bf{y}}, \pmb{\theta}, t)]$; the value obtained was 3.123 × 10⁻³ which also confirms the practical identifiability. Table 4 presents the parameter correlation matrix. No linear dependence between the parameters is found which robustifies the appropriateness of the model structure to represent the dynamics of the elastic modulus in the gelation and, consequently, of hybrid crosslinking kinetics. Figure 4 shows the agreement between model predictions and the rheological data. We note that the dispersion of the measurements increases with time which corroborates the adoption of constant relative variance model for prediction error, and ultimately using the criterion.

Although the assumptions and the eventual need of extending the experiments, the models provide a good basis to understand the undergoing kinetic transformations. Besides the quality of the model in explaining the kinetic transformation, good agreement with empirical models generally employed to represent polymerization, particularly in gelation reactions, is obtained.

5.   Conclusions

This paper introduces the application of dynamic optimization to fit the kinetics of self-crosslinking and hybrid crosslinking of chitosan-based systems in gelation process. The model describes the kinetics of the transition of the solution state to gel state, and we use rheological data (the elastic modulus) to fit empirical kinetic models to describe the dynamics of the rheological degree of conversion. The model representing the kinetics of conversion of solution to gel state has the form of a DAEs system and we adopt the Maximum Likelihood criterion for model fitting as it maximizes the freedom regarding to the parametrization of the error measurement variance model which varies over the rheological monitoring time window.

The model fitting problem is numerically solved using a sequential approach where the procedure iterates between the integration of the DAEs and the solution of the optimization problem to maximize the log-likelihood until convergence. During the integration step, the parametric sensitivity matrix is also constructed, and subsequently used to measure the model adequacy and model accuracy. Specifically, in this post-analysis step we compute ⅰ. the parameters confidence intervals; ⅱ. the parameters correlation matrix; ⅲ. assess the practical model identifiability using the FIM; and ⅳ. assess the model adequacy to explain the data using the $\chi^2$ goodness of fit test.

The numerical approach was successfully applied to fit the kinetics of ⅰ. the self-crosslinking of chitosan; and ⅱ. the hybrid crosslinking of chitosan with genipin 0.15 %(w/w). For both cases the results (i.e., kinetic rates and reaction orders) are in good agreement with results available in the literature; for a comparative analysis, see the discussion in Sections 4.1 and 4.2. The models identified are statistically adequate to represent the process data features and practically identifiable despite of their empirical nature. Further, the models are of reduced complexity and the parametric confidence intervals are tighten enough to prospectively allow attaining good performance in prediction. The reaction order obtained for self-crosslinking kinetics is 1.3375 $\pm$ (0.0151) and for hybrid crosslinking with chitosan is 2.2402 $\pm$ (0.0373). Further, for self-crosslinking $k_1$ = 5.898 × 10⁻⁵ ± 3.92 × 10⁻⁷ and $k_2$ = 8.304 × 10⁻³ ± 3.48 × 10⁻⁵, and for hybrid crosslinking we obtained $k_1$ = 8.264 × 10⁻⁴ ± 9.22 × 10⁻⁶ and $k_2$ = 4.163 × 10⁻² ± 5.04 × 10⁻⁴.

Finally, we believe that the procedure introduced can easily be applied to other systems and polymer reactions when rheological data is used for monitoring the undergoing transitions. A topic that worth being analyzed in the future is the extension of the methodology to first principles models describing the concentrations of polymer chains of different sizes or equivalent moment-based models.

Conflict of interest

The authors declare there is no conflict of interest.