Complex networks for data-driven medicine: the case of Class III dentoskeletal disharmony

A group of 70 Class III patients (Caucasian ancestry, 7–13 years of age; mean 9.5 years; 40 females and 30 males) from the Graduate Orthodontic Program at the University of Michigan, Ann Arbor, Michigan, and from the Department of Orthodontics, University of Florence, Italy, were analyzed before orthopedic treatment (T1) with rapid maxillary expansion combined with a facial mask (RME/FM). The same patients were re-evaluated (T2, 11–18 years of age, mean 14.7 years) with a lateral cephalogram (see figure 2) following a second phase of treatments with fixed appliances (braces). The patients were matched with a cohort of 70 untreated Class III controls (U2 11–18 years of age, mean 14.5 years) on the basis of race, sex, age, and type of malocclusion. The control patients were selected randomly from a cohort of 1070 Caucasian Class III subjects from the same data sources. All subjects were enrolled previously in large descriptive estimates of craniofacial growth in Class III malocclusion [33, 34] and in two clinical studies [35, 36].

Zoom In Zoom Out Reset image size

To be included in this study, the subjects had to satisfy all of the following criteria:

• Caucasian ancestry;
• No orthopedic/orthodontic treatment prior to the initial cephalogram;
• Diagnosis of Class III malocclusion based on anterior crossbite, accentuated mesial step relationships of the primary second molars, permanent first molar relationship of at least one half cusp Class III, a negative Wits appraisal ($\lt$−2 mm), and ANB angle less than 0;
• No congenitally missing or extracted teeth; and,—No craniofacial syndromes.

The three components of the RME/FM therapy used in this study were a maxillary expansion appliance, a facemask, and heavy elastics [35, 36]. Treatment was initiated with the placement of a bonded or banded maxillary expander to which were attached vestibular hooks extending in a superior and anterior direction.

Patients were instructed to activate the expander 1–2 times per day until the desired transverse width was achieved. Patients were given facemasks with pads fitted to the chin and forehead for support either during or immediately after expansion. Elastics were attached from the soldered hooks on the expander to the support bar of the facemask in a downward and forward vector, producing orthopedic force levels up to 350–600 grams per side. Patients were instructed to wear the facemask for a minimum of 14 h per day. All patients were treated at least to a positive dental overjet before discontinuing treatment; most patients were overcorrected toward a Class II occlusal relationship. The average duration of RME/FM treatment was 1.1 years.

The majority of the patients (about 90% of the cases) underwent a second phase of treatment with fixed orthodontic appliances. On average, fixed appliance therapy lasted 18 months.

A cephalometric analysis comprising N = 21 variables (10 linear and 11 angular: see table 1 and figure 2) was performed. Data from landmarks contained in each cephalogram were entered into cephalometric software (Dentofacial Planner Plus^TM, Version 2.5, Toronto, Ontario, Canada). A standardized enlargement factor (8%) was applied to all linear cephalometric measurements. The error of the method for the cephalometric measurements was evaluated by repeating the measures in 30 randomly selected cephalograms. Error was on average 0.8° for angular measures and 0.9 mm for linear measures.

The growing craniofacial system can be modeled as an aggregate structure of a variety of agents in which the clinical (e.g., radiographic, functional) characteristics can be represented as nodes of a network where the relationships between these characteristics are represented as link among such nodes. Each component of the network can be regarded as a processing unit of information. Network analysis can deepen the entanglement of the elementary components or interactions [37]. The specific entities interacting with each other depend on a type of underlying reticular structure that may have a large effect on the behavior of the system. A link can represent any sort of interaction between the various parts of a system; since we resort on angular and linear data (table 1), we can only study relations among the parts and not causality. Hence, we will take as a measure of interaction among the N = 21 cephalometric variables ${{x}_{i}}\;,\;i=1\ldots N$ the Pearson correlations between two variables:

$$\begin{eqnarray*}&&{{\rho }_{ij}}=\frac{\left\langle \left( {{x}_{i}}-\left\langle {{x}_{i}} \right\rangle \right)\left( {{x}_{j}}-\langle {{x}_{j}}\rangle \right) \right\rangle }{\sqrt{\left\langle {{\left( {{x}_{i}} \right)}^{2}} \right\rangle -{{\left\langle {{x}_{i}} \right\rangle }^{2}}}\sqrt{\langle {{\left( {{x}_{j}} \right)}^{2}}\rangle -{{\left\langle {{x}_{j}} \right\rangle }^{2}}}},\end{eqnarray*}$$

where $\left\langle \;\; \right\rangle$ indicate the ensemble average. Since we have a series of N_p (where N_p is the number of patients) measurements for each of the x_i variables, we will estimate the population Pearson correlation ${{\rho }_{ij}}$ by the sample correlation coefficient r_ij:

$$\begin{eqnarray*}&&{{r}_{ij}}=\frac{\mathop{\sum }\limits_{p=1}^{{{N}_{p}}}\left( x_{i}^{p}-\overline{{{x}_{i}}} \right)\left( x_{j}^{p}-\overline{{{x}_{j}}} \right)}{\sqrt{\mathop{\sum }\limits_{p=1}^{{{N}_{p}}}{{\left( x_{i}^{p}-\overline{{{x}_{i}}} \right)}^{2}}}\sqrt{\mathop{\sum }\limits_{p=1}^{{{N}_{p}}}{{\left( x_{j}^{p}-\overline{{{x}_{j}}} \right)}^{2}}}},\end{eqnarray*}$$

where x_i^p is the value of the ith variable observed in the pth patient and $\overline{{{x}_{i}}}=\sum _{p=1}^{{{N}_{p}}}x_{i}^{p}/{{N}_{p}}$ is the arithmetic mean of the ith variable. The elements ${{r}_{ij}}\in \left[ -1,1 \right]$ define the correlation matrix $R=[{{r}_{ij}}]$; since we have different groups of patients, we will indicate with R^T1 the correlation matrix corresponding to the T1 group (patients before the treatment), with R^T2 the one corresponding to T2 (treated patients) and with R^U2 the one corresponding to U2 (the control group of untreated patients).

The possible relations among the variables correspond to a highly-dimensional systems of $N\left( N-1 \right)/2=210$ relations; to reduce the dimensionality of the system, we filter the correlation matrix R by associating a weighted network ${{G}_{t}}\left( R \right)=(V,E,w)$ where the nodes V correspond to the N cephalometric variables, the edge set is $E=\left\{ \left( i,j \right):|{{R}_{ij}}|\gt t \right\}$ and the weight w associate to each link $\left( i,j \right)$ is the correlation R_ij among the nodes i and j. In such a way the graph ${{G}_{t}}\left( R \right)$ represents a relation among nodes only with correlation higher than t; for simplicity, we will indicate as $G_{t}^{T1}$, $G_{t}^{T2}$ and $G_{t}^{U2}$ the networks corresponding to R^T1, R^T2 and R^U2.

Many networks of scientific interest, including metabolic, regulatory, and computer networks, are found to divide naturally into 'modules' or 'communities', i.e. groups of densely associated components connected to each other with loose links [38]; weʼll see that this is the case also for orthodontic networks. In an orthodontic network, modules are sets of traits that are integrated internally by interactions among traits, but are relatively independent from other modules. More formally, given a graph $G=(V,E)$, a community finding algorithm finds a partition $P=\left\{ {{C}_{1}},\ldots ,{{C}_{k}} \right\}$, ${{\bigcup }_{i}}{{C}_{i}}=V$, ${{C}_{i}}\bigcap {{C}_{j}}=\varnothing$ of its vertices that can be considered a good community structure when the proportion of edges inside the C_iʼs (internal edges) is high compared to the proportion of edges between them.

Breaking complex networks into modules helps to simplify their structure and gives valuable insights into network organization and function. These highly interconnected sub-graphs represent local regions of optimization, auto-reinforcing structures that act as 'attractors' during the growth process of several biological organisms [39]. These sub-networks often operate as critical components of the whole network; conveying useful information, they provide evidence for a modular view of the networkʼs dynamics [40, 41]. The network decomposition into modules often is helpful for understanding the functional organization of the system; many biological complex networks are built up from several interlinked sub-graphs.

To detect the modules, we apply the walktrap algorithm [42] that highlights the dense subgraphs of sparse graphs using a measure of similarities between vertices based on random walks: loosely speaking, the walktrap splits the graph in regions where random walks tend to be trapped. Walktrap allows to find modules without having to fix the number of modules in advance; moreover, being based on random walks, it is expected to give results related to the spectral partitioning of a graph [43]. While many alternative module-finding algorithms are possible, walktrap is often faster and simpler to implement [42].

Notice that other metrics could be used in analysing the network structure of orthodontics data, like centralities and degree distribution [32] leading eventually to highlight the nodes that can be the driver variables for the growth of the orofacial system [31]; however, the strong correlation found among the modules' nodes hints that a hollistic approach is more indicated in tackling the understanding of orofacial growth.

By varying the threshold t, the networks $G_{t}^{T1}$, $G_{t}^{T2}$ and $G_{t}^{U2}$ will range from almost complete graphs for t = 0 (since any non-zero correlation is considered, all the links are present), to fully disconnected graphs for t = 1 (no links are present since on a finite sample the probability of getting a 100% correlation is highly unprobable). In this first phase of network analysis, we have resorted on visual inspection by expert orthodontists to discriminate the bordering situation where networks are simple enough to be examined by eye yet complex enough to convey meaningful clinical information. We find that in a region of $\pm 10\%$ around the threshold t = 0.5 the resulting networks satisfy such criteria on the other hand, weʼll see in section 3.2 that the modules detected by automatic community-finding algorithms confirm our empirical analysis. In all this section, weʼll show figure obtained using a threshold t = 0.5; for the sake of clarity, disconnected nodes (i.e. variables whose correlations are all $\lt 0.5$) will not be represented.

Figure 3 illustrates the correlation network of the cephalometric characteristics of the Class III patients (age 7–13 years) before orthopedic treatment with RME/FM (T1). We see that the network splits into two components: the dentoalveolar adaptive variables (figure 3, top) and the skeletal variables (figure 3, bottom). Moreover, the skeletal horizontal variables (figure 3, bottom right) and skeletal vertical variables (figure 3, bottom left) are loosely connected by a link among the gonial angle (ArGoMe) and the distance between Gonion and Pogonion points (Go-Pg).

Zoom In Zoom Out Reset image size

The same patients were re-assessed at the end of the RME/FM therapy, followed by a second phase of treatment with fixed appliances (T2; figure 4). Network analysis show a more interlinked structure between skeletal and dentoalveolar adaptive nodes and the appearance of a strong correlation among the overbite and overjet variables.

Zoom In Zoom Out Reset image size

Finally, figure 5 illustrates the correlation network of the cephalometric characteristics of the cohort of the untreated Class III controls of 11–18 years of age (U2).

Zoom In Zoom Out Reset image size

While the network of treated patients at the end of RME/FM therapy showed a highly interconnected structure between dentoalveolar adaptive and skeletal variables (figure 4), the network of control patients exhibited two separate sub-networks pertaining to skeletal (left) and dentoalveolar (right) adaptive nodes (figure 5). The core of the network of treated patients consisted of a triangle composed by three maxillomandibular divergence variables: SN-GoGN, ArGoMe, and PP-MP. These variables seem to act as a bridge between skeletal features of Class III malocclusion and the dentoalveolar adaptive components of the craniofacial system.

We detect modules by applying the walktrap algorithm [42] on the complete graphs corresponding to consider all the possible edges weighted with their correlation r_ij; in particular, we use the igraph [44] implementation of the function walktrap.community of R [45]. To compare the results with the empirical outcomes of section 3.1, we show in all the pictures only links with a correlation correlation ${{r}_{ij}}\gt 0.5$ and only nodes connected by such links. The analysis of the modules shows interesting differences among the three groups of patients: in fact, the modules characterising the networks of the patients before control (T1) and the untreated control patients (U2) share a similar structure (see figures 6 and 8) with three modules composed by analogous nodes. In fact, like in figure 3, the modules of T1 shown in figure 6 correspond to dentoalveolar adaptive nodes (top right of the picture), to skeletal horizontal nodes (bottom right of the picture) and skeletal vertical nodes (left of the picture); a similar correspondence exists among figures 5 and 8. Notice that for both T1 and U2 the modules are either separated or linked by negative links, indicating that such groups of nodes work as separate, non interacting, oro-facial modules. On the other hand, the network of treated patients (T2, figure 7) shows more modules ; apart to the new (respect to T1 and U2) isolated community formed by overjet and overbite (whose correlation is due to the successful action of the braces), the remaining communities are much more interlinked, possibly hinting a harmonization of the orofacial modules that start working together thanks to the orthodontic correction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Hence, the arrangement of modules of at the start (T1, figure 6) and at the end of treatment (T2, figure 7) of the same Class III patients confirms the strict interlinked module topology after RME/FM treatment, as compared with untreated controls (U2, figure 8).

Orthodontic practitioners evaluate dentofacial phenotypes from clinical, functional, and cephalometric data. It is likely that in the near future these large-scale, high–dimensional datasets can be integrated by a remarkable diverse framework. Ultimately, this 'organized complexity' is worth addressing when orthodontists are forced to make therapeutic choices under conditions of uncertainty, when they realize that corrective appliances will stimulate the algorithms of the living orthodontic world, discovering causality beyond therapeutic objective, unexpected consequences that seem to contradict any reasonable clinical scenario. This study showed that treatment of Class III dentoskeletal disharmony is able to induce a more interlinked network structure of cephalometric features. A more interlinked connection between adaptive and skeletal nodes could enhance the self-correction property of dentoalveolar nodes thus improving the effects of Class III treatment.

In general, our approach could be easily applied to other fields of medicine, allowing modern medical practitioners to take full advantage of the constantly increasing amount of patients' information due to the availability of novel and ever more sophisticated diagnostic procedures. In fact, network approaches to medicine not only allow to mine information both from patients' records [31, 32] and from medical literature [50], but also to understand of how organ systems dynamically interact and collectively behave as a network to produce health or disease and different physiologic states, as recently shown in [51–53].