A mechanistic model to study the effect of the cell wall on starch digestion in intact cotyledon cells

The role of the plant matrix is recognized as the main factor restricting starch digestibility in beans. Several authors have provided insights about the mechanisms behind the reduced starch digestibility in plant matrices. In this study, by means of a mathematical model, we provide a mechanistic explanation of the role played by the cell wall. It was confirmed that starch entrapped within intact cells could only be hydrolysed after α-amylase diffusion through the cell wall. This process is limited by the pores naturally present in the cell wall and the adsorption of α-amylase to the cell wall surface. These factors restrict the concentration of α-amylase available within the cells. The model assumptions are valid under controlled laboratory conditions and were validated with in-vitro digestion data giving very accurate results. The proposed approach provides new information to understand the digestibility of starch, and possibly other macronutrients, in complex food matrices.


Introduction
Starch is the primary carbohydrate produced by plants as energy storage and one of the main macronutrients in the human diet (Singh, Dartois, & Kaur, 2010). Due to the primary role of starch in nutrition, substantial research has been conducted in order to understand the factors affecting its hydrolysis during gastro-intestinal digestion, from which food structure is one of the most relevant (Parada & Aguilera, 2011).
A good example showing the importance of food structure in starch digestion are legumes, a food product with a relatively low glycaemic index. This has been confirmed by in-vitro studies, which indicate slow digestion kinetics of starch in beans (Berg, Singh, Hardacre, & Boland, 2012). Several publications confirm those results and demonstrate that the limited starch digestibility in beans is due to the restricted passage of α-amylase through the cell wall (CW) (Dhital, Bhattarai, Gorham, & Gidley, 2016;Pallares Pallares et al., 2018;Rovalino-Córdova, Fogliano, & Capuano, 2018). Even though this information provides insights about the mechanism behind the low glycaemic index of legumes, further research is needed to fully understand this effect.
In the past few years, several researchers have postulated the use of models in an attempt to describe the mechanisms that govern starch hydrolysis. For instance, Goñi, Garcia-Alonso, & Saura-Calixto (1997) modelled the digestibility of starch granules in raw and cooked conditions. Later on, Al-Rabadi, Gilbert, & Gidley (2009) and Edwards, Warren, Milligan, Butterworth, and Ellis, (2014) proposed mathematical models to estimate the level of starch digestion in grains taking into account differences in particle size. These models used only kinetic parameters, assuming that the only factors affecting starch hydrolysis were related to the enzymatic conversion of starch. This might be the case for simplified systems such as free starch granules, but when studying complex food matrices other aspects such as (enzyme) transport phenomena should also be considered. To the best of our knowledge, only one study has been proposed so far, in which a mechanistic approach has Abbreviations: CW, cell wall; ICC, intact cotyledon cells; SIF, simulated intestinal fluids; Pe, partition coefficient; Ce, enzyme concentration; Cs, substrate concentration; Ce b , enzyme concentration in bulk; Cs b , substrate concentration in bulk; Ce w , enzyme concentration in cell wall; Ce in , enzyme concentration within intact cells; Cs in , substrate concentration within intact cells; Ce av , enzyme concentration in cell wall and within intact cells; φ, volume fraction; D, diffusion coefficient; δ, cell wall thickness; Kd, hindrance coefficient; V IC , volume intact cells; V T , total volume; k, mass transfer coefficient; A, total surface area; V w , volume of cell wall; V in , volume within intact cells; k enz , kinetic constant; MDC, mechanically damaged cells. been used (Do et al., 2020). However, this study does not provide detailed information regarding the transport of enzymes through the cell wall and the factors affecting it.
Recently, it was demonstrated that α-amylase can bind to cellulose and bran fibre, thus serving as an inhibitor for amylolysis (Dhital, Gidley, & Warren, 2015). However, the real impact of such interactions in a system where starch is entrapped within a CW matrix has been rarely addressed  and deserves further investigation.
The aim of this study is to develop a mechanistic model to understand how the CW properties and its interaction with α-amylase influence starch hydrolysis in legume cells under idealized conditions. We hypothesize that the hindered diffusion of α-amylase through the pores of the CW can be explained by the synergistic effect between steric factors and adsorption interactions to CW components.

Sample preparation
Intact cotyledon cells (ICC) from red kidney beans were isolated following the procedure described by Rovalino-Córdova et al. (2018). Briefly, soaked beans were boiled for 1 h and sieved. Samples collected in a sieve with a mesh size of 70 μm were used for analyses.

In-vitro digestion experiments
In-vitro digestion was carried out using a modified version of the protocol developed by Minekus et al. (2014). In short, 2.5 g of ICC were gastric digested by the addition of 2.5 mL of simulated fluids (pH3) containing pepsin. The resulting gastric chyme was combined (1:1 v/v) with intestinal fluids (SIF), trypsin, chymotrypsin, α-amylase and incubated for 4 h at pH7. The experiment was conducted in 50 mL falcon tubes; samples were mixed by means of vertical rotation (20 rpm) and incubated at 37 • C. In this study, all experiments were performed using substantially lower concentration of α-amylase than what normally employed in conventional in-vitro digestion experiments. This was done to ensure that the substrate was always in excess so the reaction velocity will be constant (Segel, 1993). Samples were collected at different time points for determining the kinetics of starch hydrolysis and enzyme concentration. In order to stay within the range of constant reaction velocity, only a small fraction of the substrate (not higher than 5%) was consumed (Segel, 1993). Therefore, only those time points that fulfilled this requirement were considered for the analysis. Before conducting in-vitro digestion experiments, all enzyme activities were determined according to Minekus et al. (2014). The purity of α-amylase was verified by SDS-page.
In a second set of experiments, in-vitro digestion at different enzyme concentrations (0.3, 1, 5, 10 and 25 U/mL) was conducted to determine the partition coefficient between the CW and aqueous phase. Samples were incubated for 20 h to ensure equilibrium between both phases. α-amylase concentration was quantified as described in the following sections.

Quantification of starch hydrolysis
Aliquots taken from in-vitro digestion were further hydrolysed with amyloglucosidase as described by Rovalino-Córdova et al. (2018) and the corresponding glucose concentration was quantified by D-glucose assay (GOPOD method) following the manufacturer's instructions (Megazyme Inc. Bray, Ireland).

α-amylase concentration
Enzyme concentration (in the bulk and ICC) was quantified after invitro digestion using a modified version of alpha amylase assay procedure (Ceralpha method) Megazyme, Inc. (Bray, Ireland). After incubation, tubes were opened and the supernatant was separated from ICC using a Falcon TM cell strainer (mesh size 70 μm). Cells collected were weighed and re-suspended in SIF in a ratio 1:7 respectively. This suspension was left overnight under stirring conditions to break down the cells (Rovalino-Córdova et al., 2018). Enzyme concentration was determined by combining 0.1 mL of sample with Amylase HR reagent and incubated at 40 • C for 20 min. Subsequently, the reaction was stopped by the addition of 1.5 mL of stopping reagent (20% tri-sodium phosphate solution, pH 11). The absorbance was read at 400 nm against distilled water. A calibration curve was constructed by plotting alpha amylase concentration versus absorbance. This curve was linear over the concentration range 0-25 μg/mL. A linear regression was fitted for quantitation. All time points tested were performed in duplicates.

Mathematical modelling approach
A mechanistic dynamic model based on Fick's law (Taylor & Krishna, 1993) was developed for enzyme diffusion within ICC, and subsequently for the diffusion of starch hydrolysis products from the inner part of ICC towards the exterior liquid phase (bulk). Considering the high affinity of the enzyme to cell wall components (Dhital et al., 2015), the measurement of α-amylase partition coefficient (P e ) allowed us to incorporate the effect of enzyme adsorption in the model. Moreover, besides the resistance exerted by the CW, the effect of the stagnant layer in proximity to the CW was included in the model by using film theory (Taylor & Krishna, 1993). Fig. 1 shows a schematic view of the concentration profiles for the enzyme (Ce) and starch hydrolysis products (Cs).
Since the enzyme concentration in the wall (Ce w ) and inside ICC (Ce in ) could not be measured separately, an average of them is considered in the model (Ce av ), which represents the measured concentration in the pellet of intact cells. Eq. 1 shows this relation in which φ stands for volume fraction (φ w + φ in = 1).
Considering that, by definition, the partition coefficient P e is the ratio of the concentration of a compound in two immiscible phases (aqueous and cell wall) at equilibrium, we can state: By combining Eq.1 and Eq. 2 and considering that Ce w(t) ≈ P e Ce in(t) , we obtain the following relation which expresses Ce in as a function of Ce av at any time: The dynamic model is represented by a system of four differential equations (Eq. 4-Eq. 7), which together with Eq. 3 are sufficient to represent the entire dynamic system. These equations represent the change in the concentration of the enzyme and hydrolysis products (Cs) over time. Eq. 4, Eq.5 and Eq. 7 have a comparable structure, in which two terms can be distinguished at the right hand side: the resistance (first bracket) and the driving force for diffusion (second bracket). The former is the sum of the stagnant layer resistance, represented by 1/k, and the resistance of the CW (1⁄ (D/δ P e kd)), which is a more complex term that includes the diffusion coefficient (D), the thickness of the cell wall (δ), the partition coefficient (P e ) and the hindrance coefficient (kd).
The driving force, on the other hand, is the difference of the concentration in the bulk and inside ICC. Evidently, equilibrium is obtained once these two concentrations reach equal values. In Eq. 4-7, A stands for the total area of the surface of ICC, while V IC and V T are the volume of ICC and the liquid phase respectively. On the other hand, Eq. 6 contains the aforementioned structure with one additional term, which represents the production rate of the hydrolysates. A detailed explanation about the derivation of this set of differential equations could be found in the supplementary material.
The model assumes that all ICCs have the same size and spherical shape, uniform CW thickness and pore size. In order to simplify the kinetics of the enzymatic reaction, all experiments were performed under excess substrate concentration allowing us to assume that constant velocity conditions were met. It is also assumed that local equilibrium takes place at the interface between water and the CW. The parameters used for the resolution of the model are presented in Table 1.
2.2.4.1. Mass transfer coefficient (k). k was calculated using the empirical expression obtained by Brian and Hales for freely moving particles with low Reynolds number (Brian & Hales, 1969).
In which d IC corresponds to the diameter of ICC, while η l and ρ l are the viscosity and density of the surrounding liquid (water) at the experimental temperature. A detailed description of this calculation could be found in the supplementary material.

Partition (P s ) and hindrance coefficient
Products of starch hydrolysis are neutral dextrins, for which the only mechanism affecting its diffusion through the cell wall is steric exclusion (Aguirre Montesdeoca, Bakker, Boom, Janssen, & Van der Padt, 2019). In this case, P s and Kd s were estimated by approximating these molecules to a spherical geometry and by considering the cell wall pores to be cylindrical. Therefore, P s and Kd s were calculated using the equations of Dechadilok and Deen (2006): The radius of DP3 dextrin and its diffusion coefficient were previously calculated by using Stokes-Einstein and Sano and Yamamoto (1993) equations respectively. A detailed description of this calculation could be found in the supplementary material.

Kinetic constant (k enz ) determination.
The reaction kinetics are represented in our model by k enz . This parameter was calculated upon starch hydrolysis of ICC whose integrity was previously disrupted to determine the catalytic action of enzymes without having the structural constraints exerted by the CW, following the protocol described by Rovalino-Córdova et al. (2018). Structural disruption of the cells was confirmed by visual inspection under a light microscope. Hereinafter these samples will be referred as mechanically damaged cells (MDC). MDC were enzymatically hydrolyzed at different α-amylase con- Fig. 1. Schematic representation of α-amylase diffusion through the stagnant layer and CW, and the exit of the hydrolysis products towards the bulk phase. Dotted lines represent the concentration gradient through the stagnant layer and CW. The dashed red line indicates the stagnant layer. Ce b = enzyme concentration in bulk; Cs b = substrate concentration in bulk; Ce w = enzyme concentration in CW; Ce in = enzyme concentration within ICC; Cs in = substrate concentration within ICC; Ce av = enzyme concentration in the combined CW + ICC compartment. centrations (0.03U/mL, 0.3U/mL, 1U/mL). These concentrations were at least 2 orders of magnitude lower than what normally used when conducting in-vitro digestion experiments in order to ensure an excess of substrate as described in section 2.2.1. Starch quantification was performed (section 2.2.2) and only those time points, in which less than 5% of starch was hydrolyzed, were considered for k enz determination. k enz was calculated from Eq. 6 considering that in this case there were no structural barriers affecting the enzyme kinetics, therefore: The results obtained for k enz at different α-amylase concentrations are detailed in supplementary material.

α-amylase partition coefficient.
For determining the partition coefficient an adsorption isotherm of α-amylase binding to ICC was constructed using enzyme concentration data at equilibrium conditions. The calculation was done assuming that at equilibrium conditions both the bulk and the inner part of ICC have the same enzyme concentration (Ce b(eq) = Ce in(eq) ) as schematically represented in Fig. 2. Moreover, for this calculation it was also assumed that the concentration of enzyme and hydrolysis products were homogeneous within ICC (no internal gradients). By definition, the partition coefficient could be calculated upon the ratio of the enzyme concentration in two immiscible phases (Eq.2), therefore, it was necessary to first determine the concentration of the enzyme in the CW.
Since Ce av is the added contribution of the enzyme concentration in the CW and within ICC, a mass balance was proposed to determine the concentration of α-amylase in the CW: Where V w , V in and V T represent the volume of CW, within ICC and total volume (V w + V in ) respectively. The only unknown from Eq. 14 is Ce w since it is assumed that at equilibrium conditions Ce b(eq) = Ce in(eq) . V T and V in were calculated using an ICC diameter of 100 μm (Rovalino-Córdova et al., 2018) and a CW thickness of 2 μm (McEwen, Dronzek, & Bushuk, 1974), while Ce av was determined experimentally. The volume fraction (∅ w ) in Eq.1 and 3 was calculated upon the volume of a sphere considering an internal (r 3 in ) and total radii (r 3 T ) as depicted in Fig. 2.

Model validation.
Starch hydrolysis products that diffused through the CW were quantified experimentally and used to validate the model at different enzyme concentrations (0.3, 25 and 50 U/mL). For this, starch hydrolysis was expressed in terms of maltotriose concentration and the only parameter considered to influence diffusion through the CW pores was its size (0.55 nm radius). Due to the lack of electric charges in the oligomers formed upon starch hydrolysis, electrical interactions with the CW components were ruled out.

Computational
Analysis. MATLAB R2017b was used for all calculations. Parameter estimations were performed using the function 'fminsearch', and the system of differential equations was solved using the function 'ode45'.

Results and discussion
In this study, we propose a mathematical model to mechanistically describe the low digestibility of starch in ICC. Previous studies have determined that starch confined within ICC could be hydrolysed by α-amylase despite of the cells' structural intactness (Berg et al., 2012;Pallares Pallares et al., 2018;Rovalino-Córdova et al., 2018). This implies that the enzyme needs to diffuse through the CW before getting in contact with the starch granules. As consequence, enzymatic hydrolysis takes place within ICC and the products of starch digestion diffuse out of the cell towards the bulk phase after hydrolysis. We consider that the complex mechanism behind starch hydrolysis in ICC is due to the simultaneous action of several factors that will be discussed in detail in this section.

Enzyme transport within ICC
Equations 1 and 2 describe the change in enzyme concentration of ICC (Ce av ) and the bulk (Ce b ), where the former increases and the later decreases as digestion proceeds. The system of differential equations was solved by using several parameters found in literature and reported in Table 1. However, other parameters like P e had to be determined experimentally in order to describe accurately the conditions present in the system. As described in section 2.2.1, different α-amylase concentrations were utilized to obtain the partition coefficient. ICC was extensively digested to provide the enzyme enough time (~20 h) to diffuse within ICC until reaching equilibrium concentrations between the supernatant and pellet prior to quantification and construction of the adsorption isotherm (Fig. 3). As seen in Fig. 3 a linear behaviour in enzyme concentration between the bulk and ICC was found independently from the initial amylase concentration used. This indicates a direct relation between α-amylase concentration in the bulk and its affinity with the CW. Due to the low concentration of enzyme utilized in this study, a saturation of CW binding sites was not achieved which resulted in a linear behaviour of the isotherm. Under these conditions, the value estimated of P e was found to be 12.19. We are aware that at higher enzyme concentrations the adsorption isotherm will eventually reach a plateau due to saturation of the binding sites present in the CW. The analysis of those concentrations was not reported since they are out of the scope of this study.
The enzyme hindrance coefficient kd e represents the resistance of mass transfer through the CW and normally occurs due to the interactions between the diffusing molecule and the pores (Ford & Glandt, 1995). In this study, kd e was estimated by fitting the model to in-vitro digestion data that quantified the enzyme concentration in the bulk phase over time (Fig. 4). A value of 1.23 × 10 -5 ± 4.94e -6 (confidence interval 95%) was obtained using data at three different enzyme concentrations. Our estimated kd e is considered very low compared to the results of other authors who investigated the diffusion of proteins through other food matrices. That is the case of the study conducted by Fardet, Hoebler, Djelveh, & Barry (1998) on the diffusion of BSA through the protein network of pasta. They found a two-fold reduction in the diffusion of BSA in comparison to its behaviour in pure aqueous solvent. Kd e value obtained by our model decreases the diffusion of α-alpha amylase by five orders of magnitude showing a great difference to what the aforementioned authors have found. We believe that the cause of such discrepancy is the structural difference between pasta and ICC. In pasta, the gluten network has a porosity ranging from 0.3-30 μm while the porosity of the CW in ICC ranges from 3.5-5.5 nm (Brett & Waldron, 1996). Considering that the CW porosity allows the passage of proteins of around 70KDa (Li, Dhital, Gidley, & Gilbert, 2019), it is reasonable to have a higher hindrance since this value generally becomes significant when the diffusing molecule and the pore size are of comparable size (Ford & Glandt, 1995).
The reduction in the transport through constricted pores is caused by two phenomena: hydrodynamic effects and equilibrium partitioning. The former refers to the frictional drag of the molecules while the later relates to steric effects and electrostatic interactions (Robertson & Zydney, 1990). These two phenomena also explain the high hindrance that our model calculated for the diffusion of amylase through CW pores. In first place, both enzyme and CW pore size have dimensions of the same order of magnitude already creating important limitations for free (and fast) diffusion . In addition to this, the interactions between the CW and α-amylase limits further its diffusion. Dhital et al. (2015) reported the inhibition of α-amylase by cellulose and bran fibre, demonstrating a strong and rapid binding interaction between these two components. This was further confirmed by a recent publication of  who showed that the binding affinity of amylase to CW made impossible the estimation of the enzyme diffusion rate by FRAP technique. Such interaction was also present in our study and could be clearly identified when plotting the concentration of enzyme present in the bulk and ICC during in-vitro digestion. Fig. 4 illustrates that after 1 h of digestion, the concentration of α-amylase in ICC exceeded the one found in the bulk, a clear indication that an adsorption mechanism was taking place. If no adsorption phenomena occurred, α-amylase concentration in both phases would have been comparable when reaching equilibrium. This adsorption process occurred in all samples tested, independently of the enzyme concentration employed.
The partition and hindrance coefficient reflect the combined effect of enzyme diffusion through constricted pores in ICC and the adsorption of α-amylase to CW components. The hindrance to diffusion due to pore size limitations are attributable to a combination of particle-wall hydrodynamic interactions and steric restrictions. The former depends on the particle proximity to the CW so any force that influences its position affects them. Even when considering those forces negligible, the finite size of the solute restricts its access to the region near the CW affecting its flux (Dechadilok & Deen, 2006).
From Fig. 4 it is evident that adsorption also has a considerable impact in the amount of enzyme that diffuses inside ICC. This phenomena is dependent on enzyme concentration since it has been found that cellulose has a limited number of binding sites for the enzyme to attach (Dhital et al., 2015). Therefore, at very low enzyme concentrations, the CW material will deplete the enzyme from the solution causing a larger effect in starch hydrolysis. This might be the most plausible explanation of the discrepancy in the degree of starch hydrolysis reported by different research groups when studying starch digestibility in ICC.
Furthermore, in Eq.4 -Eq.7 the effect of the stagnant layer is also taken into account to describe the system. This layer is an additional diffusion barrier immediately adjacent to the CW, next to a region of slow laminar flow in which convection due to stirring do not cause any significant mixing of the solution and only diffusion takes place (Barry & Diamond, 1984). In our study, mixing conditions were applied to all treatments trying to resemble what is normally occurring during digestion. We found that the resistance opposed by the stagnant layer (ke -1 ) was 2000 times smaller compared to the resistance of the CW ((D/δ P e kd) − 1 ). Therefore, we can state that the stagnant layer has a very limited effect in delaying the diffusion of α-amylase through the CW. The impact of this layer could be more relevant in other systems with higher viscosities or reduced mixing conditions (Dhital, Dolan, Stokes, & Gidley, 2014).

Starch hydrolysis within ICC
As described before, α-amylase needs to overcome adsorption and diffusion restrictions in order to hydrolyse starch. Due to the nature of those interactions, the average concentration (Ce av ) of α-amylase measured experimentally in ICC does not reflect the amount of enzymes that were actively involved in starch hydrolysis. Based on the mathematical model, we were able to determine and quantify the amount of enzymes that were capable of hydrolysing starch (Fig. 5). In the first 30 min of digestion, the amount of catalytically active α-amylase represented 25% of the amount originally present in the bulk. As enzyme concentration reached a plateau, this value increased up to 3.98 μmol/ m 3 , that is, 75% of what present in the bulk.  In the present model, enzyme kinetics was calculated when less than 5% of the total starch was consumed (constant velocity). In this situation, the excess of substrate allows all the active sites of the enzyme to be complexed (Dona, Pages, Gilbert, & Kuchel, 2010). We determined that at an enzyme concentration of 0.3U/mL, the system remained at constant velocity for 200 min allowing us to model starch hydrolysis for a considerable amount of time. We are aware that these enzymatic concentrations do not resemble the physiological environment; however, we believe that at higher enzyme concentrations the same transport mechanism will take place.
In order to represent starch hydrolysis for our experimental conditions, it was necessary to determine the kinetic constant of the reaction (k enz ). For ICC, the presence of the cytoplasmic matrix represented a structural constraint that could influence starch kinetics. For this reason, MDC were used as substrate to reflect as consistently as possible the properties of ICC without CW entrapment. By doing so, we ensured that the enzyme could work at constant velocity but still considering the influence of the cytoplasmic matrix that was otherwise impossible to be determined independently. A k enz value of 3.28 × 10 4 min -1 was determined under this experimental conditions. We decided to include the cytoplasmic matrix effect in our kinetic constant since the presence of proteins affects starch hydrolysis due to the binding of α-amylase to insoluble proteins (Yu et al., 2018). It is important to highlight that starch in MDC was hydrolysed despite of the presence of CW material since the nature of these interactions are non-active site mediated (Dhital et al., 2015).

Diffusion of starch hydrolysis products towards the bulk phase: model validation
As starch digestion proceeds, the products formed by enzymatic hydrolysis will diffuse through the CW due to the generation of a concentration gradient between ICC and the bulk phase. From literature it is known that α-amylase produces small dextrins being maltotriose one of the most abundant (Banks, Greenwood, & Khan, 1971). This was further confirmed by HPAEC-PAD performed at different time points during in-vitro digestion (supplementary material). The hindrance coefficient estimated for maltotriose was four orders of magnitude larger than what found for amylase (0.55; 1.23 × 10 -5 respectively). This was due to the lack of electric charges in the oligomers which allowed their passage towards the bulk phase at a faster rate.
The model validation was done using data sets obtained from different experiments than those employed for kd e determination. As seen in Fig. 6, the amount of oligomers predicted by the model goes in accordance to what was determined experimentally in samples with an enzyme concentration of 0.3U/mL. At higher amylase concentrations, the model is able to predict the product formation for the time range in which the reaction velocity remains constant. After this, the depletion of substrate affects the kinetics of the system (after 20 and 10 minutes for concentrations of 25 and 50 U/mL, respectively) where the model overestimates the product formation.
A sensitivity analysis was performed in order to quantify the effect of variations in ICC diameter and CW thickness on the model prediction. For this, the model was ran maintaining all parameters constant except for the two aforementioned ones where a variation in ±30% of its size was utilized. This percentage was selected based on the size distribution of ICC found in the study conducted by Do et al., 2020. The results of this analysis could be found in the supplementary section ( Figure S6). In general, the effect of these variations was small on the model predictions, which proves the robustness of the model.
Overall, we believe that our model provides a strong starting point to understand the mechanism by which starch is hydrolysed in beans. Several authors have described this phenomena in the past but only one study has provided a mechanistic interpretation of it (Berg et al., 2012;   6. Starch hydrolysis products in the bulk phase during in-vitro digestion of ICC at different enzyme concentrations. Continuous line represents the model prediction and red crosses the data obtained experimentally. For reference, a secondary axis has been included where the % of starch hydrolysis is depicted (black dots). All experiments were performed in duplicate. Do et al., 2020;Pallares Pallares et al., 2018;Rovalino-Córdova et al., 2018). This is the first time that a mechanistic model has been developed in which it is possible to discriminate between the different factors involved in starch digestion (e.g. CW, enzyme kinetics, stagnant layer) of legume matrices. Furthermore, our model provides quantitative evidence of the role played by the CW and it allows us to understand how its properties (porosity and composition) affect the diffusion of amylase. We believe that this model represents the first step towards the development of a universal model capable of representing starch digestion in complex matrices under non-idealized conditions.

Conclusion
The consumption of plant-based diets and whole grain foods has been increasing over the last years. Therefore, understanding how those complex matrices are digested should be of great concern for food scientists and industry. This study proposed a mechanistic approach designed to investigate the physical and kinetic factors involved in the reduced starch digestibility of beans. The hindrance coefficient found in our model provides a quantitative evidence of the important role that the combined effect of the constricted pores present in the CW and its interactions with α-amylase play in the delay of starch hydrolysis. We believe that the knowledge gathered by this investigation is of significant relevance and helps to understand the complexity of the food matrix and its implications in hydrolysis. Finally, the outcome of this research has wider implications than starch digestibility in beans. Other macronutrients encapsulated in different legume sources or plant-based matrices could follow similar hydrolysis mechanisms since all digestive enzymes may be affected in a comparable way by the presence of the CW and cytoplasmic matrix.