Modeling the Diet Dynamics of Children: the Roles of Socialization and the School Environment.

Childhood obesity is a health emergency in many parts of the world including the U.S. and, consequently, identifying local, regional or national intervention models capable, of altering the dynamics of obesity at scales that make a difference remains a challenge. The fact that consumption of healthful foods among most youth has yet to meet recommended nutritional standards highlights a lack of effective policies aimed at addressing the epidemic of obesity. Mathematical models are used to evaluate the roles of socialization and school environment on the diet dynamics of children. Data suggest that standard nutrition education programs may have, at best, minimal impact in altering diet dynamics at the population-level. Inclusion of peer influence (model as contagion) reinforced by the use of culturally-sensitive school menus (environmental disruption) may prove capable of modifying obesity enhancing diet dynamics; altering the diets of a significant (critical) proportion of youngsters. A framework is introduced to explore the value of behavior-based interventions and policies that account for the sociocultural environments of at risk communities. These models capture carefully choreographed scenarios to account for the fact that when dealing with diet-dynamics systems, thinking additively is not enough as it cannot account for the power of nonlinear effects.


Introduction
The CDC estimates that 12.7 million U.S. children (or 17%) of 2-to-18-year-olds are considered obese as of 2011-2012 (Ogden & Carroll, 2010). These young individuals have higher risk of cardiovascular disease, pre-diabetes, bone and joint problems, sleep apnea, and social and psychological problems than the non-obese. It has been shown that a diet abundant in fruits and vegetables can improve health and lower the risk for obesity, cardiovascular disease, and cancer (Birt, Hendrich, & Wang, 2001;Boeing et al., 2012;Murphy et al., 2012;van Duijnhoven et al., 2009) and yet, studies have shown that 90% of children (2-to-18-year-olds) do not consume the United States Department of Agriculture's (USDA) recommended daily servings of fruits and vegetables (Center for Disease Control and Prevention, 2014).
Schools are at the frontier of diet-dynamics interventions since approximately 47% of the total daily energy consumed by children comes from participation in school lunch and breakfast programmes where vegetables are under-consumed (Byker, Farris, Marcenelle, Davis, & Serrano, 2014). Approximately 19 million schoolchildren participate in the national school lunch programme (Briefel, Wilson, & Gleason, 2009). Hence, schools provide the environment where implementing and evaluating nutrition education programmes can take place. National and state-mandated policies and guidelines such as USDA 'My Plate' do influence the types of foods available to children in 99% of public schools and 83% of public/private schools (Fox, Hamilton, & Lin, 2004). Schools, an underutilised research space, provide the ideal setting for promoting and evaluating policies and programmes that foster healthy eating behaviours. Mathematical models are introduced to shed some light on the potential impact of innovative interventions in altering diet dynamics. Children tend to eat what their friends or relatives eat and so modelling our diet as the result of a 'contagion' process among young individuals, in school settings, seems like a promising direction. Furthermore, prior studies have identified the importance of socialisation (e.g. peer influence) in school environments on eating behaviours (Ammerman, Lindquist, Lohr, & Hersey, 2002;Blanchette & Brug, 2005;Katz, 2009;Lytle & Achterberg, 1995). Here, two types of simple dynamic models are introduced as a way of assessing the role of school environments on the nonlinear dynamics of socialisation at the population level. This paper is organised as follows: a brief overview of eating behaviours in children is provided in Section 2; a mathematical framework where the role of socialisation is studied under controlled highly simplified scenarios is introduced in Section 3; the analyses of a nutrition educational programmes where recidivism is common is the focus of Section 4; the results of implementing two types of standard nutrition educational programmes ('weak' and 'strong') in the presence of recidivism are analyzed in Section 5; and most importantly the findings of a behavioural-based nutrition programme in the absence of recidivism are addressed in Section 6 and with recidivism in Section 7. Finally, Section 8 collects our conclusions and observations.

Modeling eating behaviors in school settings
Healthy eating habits must include a variety of foods, high consumption of fruits and vegetables, and low consumption of fats, that is, a balanced diet, which help in the prevention of chronic diseases. The unpleasantness of healthier foods and some culturally-acclimated food taste make the systematic consumption of healthy foods difficult, since after all eating behaviours are learned over time in social settings shaped by biological, sociocultural and psychosocial factors (Ammerman et al., 2002;Blanchette & Brug, 2005;Capaldi, 1996;Katz, 2009;Lytle & Achterberg, 1995). Health behaviour theory and models that incorporate individual characteristics (self-efficacy, genetics, taste and food preferences, beliefs, skills), physical environments (available, accessibility), social (peer/parental influence, role modelling, norms), and landscape macro-system forces (culture, policy) have been utilised to assess their role in shaping dietary patterns (McLeroy, Bibeau, Steckler, & Glanz, 1988). Our proposed models include parameters that represent individual characteristics (taste and food preferences), physical environment factors (school environment) and social environment elements (nonlinear interactions among peers). Here, we utilise social-ecological theory (see Figure 1) in building population-level models with the goal of capturing in natural ways the impact that physical-, social-, and individual-level changes, have on shaping the diet dynamics of young individuals within school environments.

Physical environment
Learning how to eat healthier foods is complex. Prior studies suggest, not surprisingly, that food intake among children is influenced by the availability and accessibility of foods in their environment (Cullen et al., 2003;Perez-Rodrigo & Aranceta, 2003;Story, Kaphingst, Robinson-O'Brien, & Glanz, 2008;Van Der Horst et al., 2007). So why, has it then been so difficult to alter the dynamics of eating behaviours? Because shifting to healthy dietary patterns requires acknowledgement of the importance of nutrition, healthy diets, a culture of eating healthy foods (systematic experiences), and continuous access to healthier foods, that is, the route to sustainable diet modifications is full of challenges and obstacles (Sallis & Glanz, 2006). In other words, lasting positive changes are tied in to our ability to alter the social and environmental landscape where food consumption takes place.

Social environment
The day-to-day interactions among peers and family members shape our eating preferences/behaviours and in the process they influence or drive our diet dynamics. At schools, the need for curricula that puts social and behavioural activities that effectively promote and increase fruit and vegetable consumption must be the norm. The use of classroom time to increase the knowledge of young people with activities that include taste-testing, cooking lessons, parental involvement, and school gardening are but some of the socialisation activities that have proved to be effective in altering, albeit most often temporarily, food preferences (Ammerman et al., 2002;Lowe, Horne, Tapper, Bowdery, & Egerton, 2004;Lytle & Achterberg, 1995). Peer role modelling has also been shown to be a significant force in altering children's eating behaviours, in other words, children's eating habits are strongly influenced by what their friends eat instead of their families, unless strict parental rules are in place (Reicks et al., 2015). Corresponding social interactions have generated increases in food association and food preference learning by providing opportunities for exposure, tasting (not just smelling or seeing) and engagement in positive social experiences during food consumption activities (Birch, 1987;Birch, Zimmerman, & Hind, 1980). Further, whenever reinforcing behaviours are modelled by parents at home, it has been observed that children's diet dynamics have experienced increases in the intake of fruits, fruit juice, and vegetables (Savage, Fisher, & Birch, 2007).

Individual factors
The 'unpalatability' of healthier foods, particularly vegetables, and the onset of neophobia, or the fear of trying something new is known to influence children's food choices, lowering dietary variety and the consumption of healthy foods such as fruits, vegetables, and meats (Cooke, 2007;Dovey, Staples, Gibson, & Halford, 2008;Howard, Mallan, Byrne, Magarey, & Daniels, 2012). Studies have shown that familiarising children with varied foods can improve the liking for and intake of novel foods, the result of repeatedly exposing them to a new food at an early age (preschool and school-aged children). Many approaches have been successful or promising when dealing with youngsters, for example, pairing a new food repeatedly with an already liked food (associative learning) can facilitate the appreciation and consumption of targeted foods (Capaldi, 1996;Wadhera, Capaldi Phillips, Wilkie, & Boggess, 2015).

Mathematical framework
Models are formulated under the premise that the social adoption of healthy and unhealthy diets can be captured/modelled, at the population-level, as a social-contagion process. This social-contagion approach has been applied, in settings where social dynamics are intense, including in mathematical studies of the dynamics of bulimia (González, Huerta-Sánchez, Ortiz-Nieves, Vázquez-Alvarez, & Kribs-Zaleta, 2003), eating behaviours (Murillo, Safan, Castillo-Chavez, Capaldi-Phillips, & Wadhera, 2016), and obesity (Ejima, Aihara, & Nishiura, 2013;Evangelista, Ortiz, Rios-Soto, & Urdapilleta, 2004;Frerichs, Araz, & Huang, 2013;González et al., 2003;Jódar, Santonja, & González-Parra, 2008) at the population-level. To our knowledge, our previous work is the first to study the diet dynamics of children in school environments (Murillo et al., 2016), thus, we expand our previous models to evaluate alternative scenarios of socialisation in schools and its impact on children's diet dynamics. Models are used to generate highly simplified controlled scenarios in order to evaluate, within the confines of this artificial set up, the impact of (school) programmes aimed at youngsters. Programmes whose aim is to improve on the quality of food choices.
A typical school population may be assumed to be composed of two types of students: 'moderately' healthy individuals, M(t), those who eat 25-50% of the USDA recommended amount of fruits, fruit juice, or vegetables (FJV) and 'less' healthy ' , L(t), that is, those who eat less than 25% of the USDA recommended amount of FJV. Since prior studies have shown that standard nutritional programmes (health education and awareness, communication skills, and skill-building) have improved diet dynamics with high levels of recidivism being the norm (Pyle et al., 2006). Communication skills and skill-building improve nutrition education in contexts where cultural relevance is important and need to be addressed so that children can acquire understanding and skills to strengthen healthy eating habits (Perez-Rodrigo & Aranceta, 2003). Further, communication skills are essential for enhancing children's competence as informed consumers to be able to make their own food choices, and understand food preparation, preservation, and storage, which are all critical skills to consume healthy and perishable foods such as FJV, which is the focus of this research (Perez-Rodrigo & Aranceta, 2003) We use this knowledge (high levels of recidivism) as the starting point for the most basic model. Next, we proceed to add and explore the roles of socialisation and the physical environment on diet dynamics within the specified modelling scenarios. In building the extended model, we make use of what we have learned on food association and food preference from the results of our pilot study (Murillo et al., 2016). We tested whether or not conditioned food-preference learning, when incorporated into a school-based intervention programme, can increase fruit and vegetable consumption in school environments. So, first, we build a model that accounts for the diet dynamics of children under standard nutrition educational programmes (Models 1 and 2 in Sections 4 and 5, respectively) and calibrate their qualitative behaviour taking into account that they are passively implemented and highly correlated with high levels of recidivism. These models then incorporate behavioural-based programmes (Models 3 and 4 in Sections 6 and 7, respectively). The expanded models account for the role of positive peer influence and/or food association and food preference learning. The models are capable of generating long-lasting healthier eating behaviours (see Table 1) but do not eradicate obesity.

Standard education setting with temporary recovery and recidivism
A simplified model that involves two eating behavioural classes is introduced to illustrate a modelling philosophy that sees the role of transitions between behavioural classes as a consequence of interactions between individuals (contagion) in fixed environments (see (Murillo et al., 2016)). The population of youngsters (of fixed size N) is divided in two classes 'moderately' healthy (M) and 'less' healthy (L) all enrolled in schools. It is assumed that a typical Pre-Kindergarten to 8th grade school student spends roughly 10 years in the same school (1/μ = 10 years). It is also assumed that standard nutrition education programmes are implemented in this setting and that they are effective at improving food choices at the per capita rate (φ). It is further assumed that L-eaters may become M-eaters as a result of the programmes defined by φ. M-eaters can shift back to L-eaters, due to recidivism. The contagion power of L individuals would be considered 'successful' at worsening food choices as long as the interactions between M and L lead to an increase in the number of L's. A flow chart diagram of the M and L, eating dynamics in this deliberately standard and simple model, is shown in Figure 2. The description and values of model parameters are collected in Table 2. Specifically, the M−L system is given by the following system of nonlinear equations, Figure 2. Standard education with temporary recovery and recidivism. Less healthy eaters, L and moderately healthy eaters, M, make up the total population. λ = βL/N, where β represents the per capita peer influence rate, φ the exposure rate to nutrition programmes, and μ the per capita entry/removal rate. where λ ≡ βL/N, models M-eaters frequency dependent interactions with L-eaters. Here, β denotes the success rate of worsening eating habits of M-eaters with λ denoting the total transition rate of M-into L-eaters. The population is not closed and so the average entry rate per school year is denoted by = μN; chosen so that the total population attending school remains constant and equal to N. The rescaled system of equations for the variables X = M/N and Y = L/N, reduces to where X+Y = 1. Parameters used in Model 1 were estimated based on preliminary data from our pilot study (see Murillo et al., 2016). The initial values for M and L were chosen based on baseline measures, which indicated that at least 90% of children ate fruits and vegetables (M-eaters), and only 10% (L-eaters) did not. Study observations and interviews with school cafeteria staff revealed that students select and consume vegetables at least 1 to 2 times per week out of the 5 days spent in school, due to the required FJV selections in the school cafeteria and under the present traditional nutrition education. This implies that the average rate at which students consume FJV per year is φ = 0.5. Therefore, taking φ = 0.5, we varied β in simulation studies to determine at which values do L-eaters approach 70% in the long-term, which we measured to be the maximum L-eater population at the end of our pilot study. We found that that β > 1.76 led to the steady-state that matched our pilot study data. The final model parameters used to produce Figure 3 are shown in Table 2. The rescaled model supports the diet problem-free state and the diet endemic state, where R c,1 (φ) = β/(φ + μ) > 1, denotes the control reproduction number, a function of φ, under which the equilibrium E 1,1 is biologically meaningful. R c,1 (φ) gives the number of secondary conversions from M to L, over 1/(μ + φ) years. From Equation (2), we obtain and since 1/μ represents the average time spent in the school and 1/φ is the average time spent in the L-class, so 1/(μ + φ) represents the total average time spent in the school as an L-eater; that is, the period of time available to influence the L-eater's food choices with those of M-eaters. We see that as φ grows, R c,1 (φ) decreases. In this simplified model, it is assumed that M-eaters do not influence the eating habits of L-eaters. It is further assumed that an L becomes an M exclusively as a result of the intervention modelled by φ. Systematically implementing these programmes have a minimal impact (dashed and dotted lines) on the distribution of Mand Leaters in the population (e.g. φ must be large in order to reduce the L-eaters population). Parameters used can be found in Table 2. E 0,1 is globally asymptotically stable if and only if R c,1 (φ) ≤ 1 while E 1,1 is globally asymptotically stable whenever it exists (i.e. when R c,1 (φ) > 1), see Appendix I. In this setting, permanently reforming the L population, requires a reduction of the R c,1 (φ) to a value less than or equal to 1. Hence, the shorter the average time spent in the L-state, the most likely that the deficient diet problem is eliminated at the population level.
The model is not used to highlight the effectiveness or lack thereof of 'purely' educational programmes on altering the prevalence of L-eaters but rather to show the role of contagion on the persistence or elimination of the L-subpopulation. The absence of nutrition educational programmes (φ = 0) yields the basic reproduction number That is, the control reproduction number becomes the basic reproduction number which is given by, the product of the effective conversion rate per L (β) and the average time a student remains in the education system (1/μ). Reaching a sociocultural environment composed of mostly M-eaters requires that, φ be large and β be small. We see that whenever R c,1 (φ) ≤ 1, the population decreases, eventually becoming composed of only M-eaters (see Figure 3).

Weak and strong education with recidivism
The nutritional programme introduced via φ is supposed to be fixed and omnipresent. In general, nutrition education programmes vary in effectiveness. Strong educational programmes are those that generate relatively fast increases in the population of M-eaters. On the other hand, it is assumed (or defined) that weaker educational programmes promote slower changes in the size of M. Under strong and weaker intervention programmes, the population is vulnerable to recidivism, the assumption behind our basic model. Next, we add one more level of heterogeneity by incorporating a third class of individuals, H(t), those that consume in their 'healthy' diet, 50% or more of the USDA recommended FJV. In this expanded model, students are exposed to strong and weak nutrition education programmes. Students experiencing the 'weaker' curriculum (δ 0 ) can shift M-to H-eaters and a 'stronger' curriculum (δ 1 ) shift L-to H-eaters. Recidivism is incorporated at two levels: H-can shift to M-eaters (α 0 ) and H-can transition to L-eaters (α 1 ). The model assumes that all individuals 'eat' and cohabit the same environment. In other words, we do not consider dynamics in multiple schools by varying parameters, instead the baseline parameters for each model remain constant to represent a single school. The constants 1/α 0 (weak) and 1/α 1 (strong) denote the average recidivism time of M-and L-eaters, respectively. This is summarised in Table 3. The model diagram is in Figure 4 and parameters described in Table 4. The M−L−H model is governed by the following system of nonlinear equations,  Here, λ = βL/N, where β denotes the per capita peer influence rate. The efficacy of a strong or weak education programme is denoted δ 1 and δ 0 , respectively. Recidivism is a possibility following the education programmes and is represented by α 0 (weak) and α 1 (strong). The parameter μ represents the per capita entry/removal rate.
where = μN and λ = βL/N. Here, 1/δ 1 is the average length of effectiveness of the strong education programme on L; 1/δ 0 represents the average period of effectiveness of the weak education programme for M-eaters. Parameters used in Model 2 were chosen based on the particular scenario that motivated this model. For the initial conditions, we assumed that the L-eater population remained the same (e.g. 10%), and that the remaining students were either H-or M-eaters. The parameter value for β was estimated using the control reproduction number from our pilot data, where β ≈ 0.33 (see Murillo et al., 2016). Since Model 2, includes H-eaters in the population, then we assume β = 0.1, which is lower then the estimated value. Since the minimum frequency of FJV consumption in the school was 1 to 2 times per week (discussed in Model 1), we slightly increase the rate (δ 1 ) for children transitioning from L to H given that the FJV eating standards are higher for H (e.g. ≥ 50% consumption of FVJ). During our pilot study, we observed that M-eaters were less likely to eat significantly more FJV since it was already part of their diet, and we also observed that a small number of children would eat additional FJV, and thus assumed δ 0 < δ 1 . To consider the scenario when the recidivism is lower for the at-risk population due to stronger education programmes aimed at the at-risk group (L-eaters), we assumed α 1 < α 0 . The final model parameters used to produce Figure 4 are shown in Table 4. We explore three recidivism scenarios: α 1 = 0 (no recidivism from H to L); α 1 > 0 (recidivism from H to L); and α 1 = rβL/N (nonlinear per-capita recidivism from H to L). The system is rescaled as follows:
The distribution of M-, L-, and H-eaters is shown in Figure 5, when model parameters are held constant, α 1 = 0 and α 0 is varied (four values). These results show that even a smaller rate of recidivism (α 0 = 0.15) can alter the distribution of eaters within a set proportion of L-eaters. That is, the proportion of L-eaters in the population N.
One can easily evaluate the relative ratio of L-eaters to the combined L-and H-eaters, that is, Figure 5. The distribution of M-, L-, and H-eaters are shown when α 1 = 0 and α 0 (no recidivism from H to L) is varied (four values). These results show that even a smaller rate of recidivism from H to M (α 0 = 0.15) alters the distribution of eaters. Parameters used can be found in Table 4. Figure 6. The relative ratio of L-eaters to the combined Land H-eaters is shown as a function of R c,2,1 denoted R 0 here. Although the control reproduction number increases, the relative ratio reaches a threshold demonstrating that the increase in R c,2,1 (α 0 , 0) will not have significant impact on the proportion of L-eaters in the population. Parameters used can be found in Table 4.
as a function of R c,2,1 (α 0 , 0) ( Figure 6). This ratio is biologically relevant when R c,2,1 (α 0 , 0) > 1, which is equivalent to β > (μ + δ 1 )(1 + δ 0 /(μ + α 0 )). When the average total of school exit rate and 'strong' education rate is smaller than the rate of conversion of M-to L-eaters (μ + δ 1 < β), then the above condition holds. A fraction of L-eaters will persist in the population if the success of 'strong' and 'weak' education programmes are negligible. Therefore, without recidivism we have the possibility of obtaining a diet problem-free state and are able to define a control reproduction number that provides the threshold value indicating either a diet problem-free or endemic state. That is, in the nutritional context, an education programme that is so effective such that there is no recidivism could potentially yield a population with a low number of L-eaters when the threshold value R c,2,1 (α 0 , 0) < 1. Similarly, even when the education programme is moderately effective at improving eating habits, the fact that there is no recidivism in this case, we observe that the proportion of L-eaters does not increase even as R c,2,1 (α 0 , 0) increases. Next, we consider the case with α 1 recidivism.

Case 2: chance of breaking 'good' diet (α 1 > 0)
If we assume the recidivism rate (or chance of breaking 'Good' diets and returning to the Leating state) to be a positive constant, or α 1 > 0, then the model has neither a diet problemfree equilibrium (E 0,2,2 ) nor a control reproduction number (R c,2,2 (α 0 , α 1 )). However, it always has a diet problem-endemic equilibrium E 1,2,2 = (X 1,2,2 , Y 1,2,2 , Z 1,2,2 ) where Figure 7 shows the equilibrium proportion of students consuming a 'less' healthy diet, Y * = Y 1,2,2 , as a function of the successful social-interaction rate β. Model parameter values used are collected in Table 4, making use of two α 1 values (α 1 = 0.09 and α 1 = 0.00009). Figure 7 shows that even for β = 0, there is an endemic proportion of L-eaters Y * = α 1 δ 0 /(α 1 (μ + δ 0 ) + (μ + δ 1 )(μ + α 0 + δ 0 )) > 0. This initial proportion of Y * is monotonically increasing in α 1 . It approaches zero as α 1 tends to zero. Moreover, Figure 7  shows that the qualitative behaviour of Y * as a function of β is significantly affected by the value of α 1 ; higher values of α 1 generate increases in the level of Y * , while low levels of α 1 imply slow increases in Y * at the start, i.e. there is a beginning phase of β for which the proportion Y * remains close to zero and then it increases steadily.

Case 3: peer influence in breaking 'good' diet (α 1 = rβL/N)
The case α 1 = rβL/N is of interest as it may help understand the possible effect of negative peer influence in shaping 'healthy' diet patterns. Here, r denotes the relative probability of 'success' of peer-influence contacts that L-eaters have on H-eaters, with respect to those that L-eaters have on M-eaters. A successful peer-influence contact means that L-eaters successfully convert H-or M-eaters into L-eaters. Numerical simulations for selected values of α1 are shown in Figure 8. The model has the same diet problem-free equilibrium E 0,2,3 (as in the case α 1 = 0 or equivalent to E 0,2,1 ). However, the control reproduction number is It is computed as the ratio between the influence rate β and the critical peer influence rate corresponding to a zero endemic level of L-eaters β 0 = (μ + δ 1 )(μ + α 0 + δ 0 )/(μ + α 0 + rδ 0 ), (Safan et al., 2006). Similar to R c,2,1 (α 0 , 0), here the R c,2,3 (α 0 , 0) includes the product of β, which denotes the rate that M-eaters transition into L-eaters, and 1/(μ + δ 1 ), defined as the total average time spent in the school in the L-state before transitioning into H. Also, the rate of 'weak' education programme with rδ 0 and without δ 0 peer influence, over the average time spent in school as an H-eater before becoming an M-eater. Further, the model supports the existence of a diet problem-endemic equilibria E 1,2,3 = (X 1,2,3 , Y 1,2,3 , Z 1,2,3 ), where X 1,2,3 = 1 − Y 1,2,3 − Z 1,2,3 , Z 1,2,3 = 1 1 − r 1 − Y 1,2,3 − μ + α 0 + rδ 0 μ + α 0 + δ 0 · 1 R c,2,3 (α 0 , α 1 ) (9) with Y 1,2,3 being the feasible solution (positive) of the quadratic equation with Y 1,2,3 ∈ [0, 1]. Equation (10) could be seen as a bifurcation equation in the parameter R c,2,2 and the variable Y 1,2,3 . The plane (R c,2,2 , Y 1,2,3 ) has the bifurcation point P 0 = (1, 0), at which the bifurcation direction is to be investigated. On applying the same procedure analysis shown in (Safan & Dietz, 2009) and making use of the implicit function theorem, it is shown that the model supports the existence of backward bifurcation whenever δ 1 > δ 1 and r 1 < r < r 2 where , If neither of the inequalities in (11) hold, then the model has a forward bifurcation (see Figure 9), that is, a unique diet problem-endemic stable equilibrium exists if R c,2,3 (α 0 , α 1 ) > 1; while no positive diet problem-endemic state exists if R c,2,3 (α 0 , α 1 ) < 1. Further, the diet problem-free state is stable. The different cases for the existence of the diet problem-endemic states are below: (1) If at least one of the inequalities in (11) does not hold, then the model of interest shows the existence of a forward bifurcation. Hence, a unique diet problemendemic state exists if and only if R c,2,3 (α 0 , α 1 ) > 1. This equilibrium is given by  Table 4.

Food association learning with positive peer influence
Prior studies have suggested that food association learning techniques have been effective at exposing and increasing the liking for vegetables at the individual-level. Here, we explore the impact of food association learning, as a contagion process, under behavioural-based intervention programme aimed at reducing L/N. We introduce a new class of food association learners, A(t), that is, eaters who learn to pair new food with an already like food to facilitate the appreciation and consumption of targeted foods (Capaldi, 1996;Wadhera et al., 2015). For example, learning to associate celery with peanut butter to consume more celery is a form of food association learning. The population level model is built under the assumption that students do not change their physical eating environment. That is, students who remain in school do so under the same set of parameters (e.g. same rate of social interactions, exposure to food association learning programmes, etc.) and these socioenvironmental factors do not change. The model is thought to involve two types of contagion/conversion processes: 1) negative where M-eaters convert to L eaters (λ 1 ) and 2) positive where M-or L-eaters convert to A-eaters (λ 2 ). M-eaters can shift to A-eaters directly via the 'strong' effects of food association learning (γ 1 ) or through social interaction with A-eaters (λ 2 = β 2 A/N). A second exposure to a food association programme for L-eaters, γ 2 , the result of the 'positive' social influence resulting from interactions between L-and A-eaters may lead to increases on A-eaters. Although γ 1 and γ 2 can be interpreted as the reciprocal of recidivism, here, recidivism is not considered explicitly in the model. Only the impact of 'positive' peer influence is Figure 10. Food association learning with positive peer influence model. Less healthy eaters, L, moderately healthy eaters, M, and food association learners, A, are included. Here, λ 1 = β 1 L/N and λ 2 = β 2 A/N, where β 1 is the per capita peer influence rate shifting an Mto an L-eater and β 2 is the per capita peer influence rate shifting an Mor L-eater to an A-eater. Here A-eaters represent a new class of individuals, who learn to pair targeted foods (fruits and vegetables) with an already like food to facilitate the consumption of fruits and vegetables. The efficacy rate of a first or second exposure is denoted γ 1 and γ 2 , respectively. Parameter μ represent per capita entry/removal rate. Lastly, r represents the relative susceptibility of L-eaters who shift to an A-eater.  Figure 10 with parameter values and definitions are listed Table 5. This new model is given by the following system of nonlinear equations, where N = M+L+A, where λ 1 = β 1 L/N and λ 2 = β 2 A/N. We rescale this model by putting X = M/N, Y = L/N, and Z = A/N. Parameters used in Model 3 were also based on our pilot study data (see Murillo et al., 2016). The initial conditions were similar to those in Model 2. The rate shifting M-to an L-eater (β 1 ) were the same as in Model 1, and we assume the rate at which M-shift to an A-eater (β 2 ) is three times less then β 1 since 'food association' learning techniques indicate that 'successful' learning requires repeated exposures, and thus β 2 < β 1 . In our pilot study, those who were identified as an L-eater were less likely to try novel pairings of vegetables with other foods, and during the intervention and testing periods, only one-third of children ate vegetables on a given day during the week. Therefore, the rate at which individuals shifted from L-to an A-eater during the year was estimated to be γ 2 = 0.06 (or once per month). Children who were already eating their vegetables, M-eaters, were more likely to consume those foods and their increase was about one-third on a given day. Hence, the rate at which M-shift to the A-state is γ 1 = 0.35 (or at least 7 times per month or 1-2 times per week). The final model parameters used to produce Figure 10 are shown in Table 5.
The rescaled model has a diet problem-free equilibrium where
Here R c,3 represents the product of β 1 , the rate of M-eaters transitioning into L-eaters, proportion of M-and L-eaters at steady-state (1 − Z 0,3 ), and the average time L-eaters stay in the school system before they transition to A-eaters with (rβ 2 ) and without (α 2 ) peer influence.
If Y 1,3 > 0 at equilibrium, then the analysis shows the existence of a diet problemendemic equilibrium whose components are given by with Y 1,3 being the solution of the quadratic equation with It is shown below that equation (20) has a unique feasible equilibrium solution.

Uniqueness of the diet-problem endemic state
To prove the uniqueness of the diet-problem-endemic equilibrium for the rescaled model it suffices to prove that equation (20) has a unique positive solution. To this end, we rewrite equation (20) in the equivalent form with It is clear that both functions G 1 and G 2 are quadratic with positive coefficients. Both functions are parabolas. They are always increasing for Y 1,3 > 0 but with different tangents. Thus, they either intersect at one point guaranteeing a unique solution, or they never intersect, that is, they have no solution. This completes the proof.
Consequently, the rescaled model has a unique diet-problem-endemic equilibrium that exists if and only if R c,3 (γ 1 , γ 2 ) > 1. Hence, if we are to apply a control strategy based on subjecting M-eaters to a strong effect of food association, then the minimum rate at which M-eaters must be shifted to A-eaters status needs to satisfy the inequality γ 1 > γ c 1 , where

Food association and food preference learning
We evaluate the impact of food association learning on developing sustainable healthy eating habits in the long-term via food preference learning (found in Murillo et al., 2016).
Here M-eaters and L-eaters will enter the food association learning programme at rate γ 1 and γ 2 , respectively. After completing the programme, a portion (p) will become food preference learners P(t) at the rate pα, in which we consider the food association learning programme successful. That is, P-eaters are those who have developed a preference for specific repeated pairings of foods, which lead to the long-term consumption of targeted foods since the novel flavour and foods have become part of the portfolio of preferred food types (Capaldi, 1996;Wadhera et al., 2015). Figure 11. Food association and food preference learning model. Less healthy eaters, L, moderately healthy eaters, M, food association learners, A, and food preference learners, P, are included. P-eaters are individuals who have developed a long-term preference for the consumption of targeted foods. The per capita peer influence rate shifting an Mto an L-eater is denoted β, where λ = βL/N. The entry rate into a food association learning programme is denoted γ 1 for M-eaters and γ 2 for L-eaters. Given the effectiveness of the programme, α, some proportion p will become preference learners. Recidivism from the Lto M-state is possible at rate φ. Parameter μ represent per capita entry/removal rate. Recidivism of A-eaters, where they return to old ways of eating as M-eaters, occurs at the rate (1 − p)α. A-eaters can shift to L-eaters after interacting with a proportion of Leaters (L/N) at rate rλ. The M-eaters who do not enter the food association programme would either maintain current eating habits or shift to an L-eater at rate λ after interacting with L-eaters. Finally, L-eaters can shift to M-eaters at rate φ. This is shown in Figure 11 with parameters defined in Table 6, also found in (Murillo et al., 2016). This next expanded system is governed by the following equations, where N = M+L+A+P is the total population size and = μN is the per-capita school entry rate. Parameters used in Model 4 were similar to those in Model 3 with a few exceptions. The initial conditions for A and P were made under the assumption that children learn to associate (A) certain foods and have developed food preferences (P) outside of the school environment. In our pilot study, 'food preference' learning was not observed since we observed children for only 8 days (see Murillo et al., 2016). A longitudinal study with additional observations in school lunch periods and surveys would be required in order to observe 'food preference' learning. Hence, in this model, some parameters such as α and φ were hypothesised. Similarly, the model (35) is rescaled by assuming X = M/N, W = A/N, Y = L/N and Z = P/N. Analyzing the model at equilibrium shows that there is a diet problem-free equilibrium E 0,4 = (X 0,4 , W 0,4 , 0, Z 0,4 ) where This equilibrium is locally asymptotically stable (see appendix A.2) if and only if the control reproduction number R c,4 (γ 1 , γ 2 ) < 1, where Moreover, the diet-problem-endemic equilibria exists and is given by E 1,4 = (X 1,4 , W 1,4 , Y 1,4 , Z 1,4 ), where and Y 1,4 is the attainable solution of the quadratic equation where Equation (37) could be considered a bifurcation equation in the plane (β, Y 1,4 ), with β as a bifurcation parameter. It has a bifurcation point (β 0 , 0), where at which the direction of the bifurcation is to be evaluated. An application of the implicit function theorem shows that a backward bifurcation exists if and only if the following set of inequalities is held where The bifurcation diagram for this model is shown in Figure 12. The solid curve represents the bigger solution (Y 1,4,1 ) and the dotted curve represents the smaller solution (Y 1,4,2 ) in the plane (R c,4 , Y 1,4 ), which exists for values of R c,4 (γ 1 , γ 2 ) < 1. Moreover, as R c,4 (γ 1 , γ 2 ) decreases, both curves get closer to each other until reaching the turning point (Safan et al., 2006) at which both of them coalesce. The value of the control reproduction number at this turning point is given by R 1 c,4 , where and In fact, the value R c,4 (γ 1 , γ 2 ) = R 1 c,4 is a threshold value, as it separates between nonexistence and existence of diet-problem-endemic states.
If we assume now that at least one of the conditions in (38) is broken, then the model supports the existence of forward (supercritical) bifurcation in the sense that a unique diet-problem-endemic equilibrium exists and is stable for R c,4 (γ 1 , γ 2 ) > 1, while no endemic state exists for R c,4 (γ 1 , γ 2 ) < 1. Hence, R c,4 (γ 1 , γ 2 ) = 1 is the threshold level separating nonexistence and existence of diet-problem-endemic states. We summarise the above results in the following: The critical control reproduction number below which diet-problem-endemic equilibria do not exist is given by if the bifurcation is forward.

Diet-problem containment possibility
One of the most important questions is how likely is the possibility to contain (get rid of) the diet problem. Here, we discuss the existence of necessary and sufficient conditions required to eliminate the diet persistence problem, with a strategy based on the application of a food association programme whose effectiveness is p ∈ [0, 1]. In the literature of mathematical epidemiology, the basic reproduction number R 0 is a key concept and the cornerstone in determining the minimum effort required to eliminate an infection whenever the model doesn't support the existence of multiple endemic equilibria. In this case, reducing R 0 below one ensures an effective control of the infection. However, over the last two decades, models that exhibit bistable endemic states, backward bifurcations and hysteresis phenomena have been introduced to the field of mathematical epidemiology. In such cases, reducing R 0 to a value slightly below one is a necessary but not sufficient condition to eliminate the infection. However, in a prior model that supports a backward bifurcation, it has been shown (see Safan et al., 2006) that the ratio R 0 /R 0 can be interpreted as a reproduction number or ratio. Hence, reducing this ratio below one ensures an effective control of the 'infection'. Thus, if we solve the inequality R c,4 (γ 1 , γ 2 )/R c,4 < 1 in terms of the probability p, we get p > p = p 1 if the bifurcation is backward, p 2 if the bifurcation is forward (41) where Formula (41) determine the critical probability (p ) of effectiveness of a diet-problem food association programme above which the diet-problem endemic state(s) disappear.  Figure 13 shows the critical level of the food association effectiveness p as a function of the contact rate β. The vertical line β = β − separates between nonexistence and existence of backward bifurcation. Therefore, for β ≤ β − , the curve p = p 2 separates between existence and nonexistence of diet-problem endemic equilibria. Thus, a probability of effectiveness slightly above p 2 ensures an effective control of the diet-endemic problem. However, if β − < β < β + , then backward bifurcation exists and p = p 1 is the threshold above which diet-problem-endemic equilibria do not exist. Thus, a food association  Table 6. This shows that even for a small p (or 0.25 < p < 0.5), the sociocultural environment can change and is sufficient for reducing the distribution of L-eaters in the population. In other words, as the proportion of P-eaters increases (p), the distribution of M-, L-, and A-eaters in the population also change, which reflects changes in eating behaviours (e.g. the sociocultural environment) of the school population. For example, notice that L-eaters eventually diminish over time when p = 0.5. programme with probability of effectiveness slightly higher than p 1 exhibits a die-out of the diet-endemic problem, where Here, the level β = β + represents the value at which p 1 hits the upper bound p = 1. Thus, for β > β + , there is no feasible value of p that ensures a wash out of the diet-endemic problem and we may seek another control strategy to first reduce the contact rate β to below β + and then apply a food association programme with high enough probability of effectiveness. This ensures an effective control of the diet-problem. Figure 14 shows a time series analysis for fixed β and different levels of p. The solid curve corresponds to the lowest value of p = 0 and the dashed-dotted solution corresponds to the largest value of p = 1. The less healthy proportion of individuals approaches zero in case of p = 1.

Discussion
In this manuscript, we have attempted to evaluate the roles of socialisation and school environments on the diet dynamics in children's school settings. We introduced four models in order to evaluate different combinations of changes in social interactions and school environments in response to the implementation of different school nutrition education programmes and policies. In Model 1, we found that some level of standard education was better than none at all albeit not that effective (when compared with those in later models). In Model 2, we found that even a low rate of recidivism lead to a constant proportion of L-eaters, which increases under the presence of negative peer influence. Model 3 has been used to show that positive peer influence and food association learning can significantly modify the culture of eaters and thus, promote healthier eating behaviours among children. Lastly, in Model 4, we found that sustained healthy eating behaviours could be achieved by conditioned food preference learning.
The results are to some degree as expected. However, from these models we have learned a lot more. For example, that non-constant peer pressure as in Model 2 (strong education) can add uncertainty, supporting two levels of L eaters but also supporting long-lasting communities that include a proportion of M and H eaters. Furthermore, we learned the importance of initial conditions as where we end up on the distribution of L, M and H eaters depends on where we start. Hence, implementing effective programmes for individuals that immigrate into a school system may be very important. In Model 4, achieving food association status (studied with Model 3) was not enough. Further, the implementation of programmes that move individuals into a new level, preference learning, is likely to have a long-lasting positive effect. Again, the role of initial conditions has been addressed since models with multiple interventions tend to support more than one distribution of L, M, A and P, which is shown in the analysis of the stability of steady states, the control reproduction numbers, as well as the bifurcation analysis and simulation studies where parameters were varied. In short, intervention efforts that account for nonlinear dynamics generate distributions of L, M, A and P that reduce the proportion of individuals in the L class while reducing the impact of recidivism.
In short, we have learned that interventions of multiple types that have proved to be effective in field studies when implemented in simple models that account for peer pressure (model as contagion), presumably instigated by education and other forms of intervention can generate sustained change and unexpected results, like multiple steady state distributions, a situation that makes it possible to reduce the problem, the size of the L class, but creating communities where further improvement is nearly impossible (backward bifurcation). We can therefore conclude that unless students that have changed their behaviours are not 'moved' into a different setting (an environment with different cultural norms) then the problem once established will persist 'regardless' of our efforts.
These frameworks provide a means for evaluating the roles of social and environmental factors from a dynamical systems approach and shed some light on the complexity of understanding eating behaviours among children in schools. However, one of the limitations of this work is the lack of heterogeneity considered in our model assumptions. Future work would incorporate individual differences in children's and their families diet, culture, and lifestyle. Similarly, children's individual response to particular programmes and the efficacy of intervention strategies should be considered in future studies. Lastly, integrating data from field studies to better assess these hypothesis-driven scenarios could be incorporated in future studies to give more accurate predictions.