Modelos bioinformáticos y estudio de receptores de proteínas mediante el uso de redes complejas para el desarrollo y diseño de fármacos eficaces en patologías del sistema nervioso central

Modelos bioinformáticos y estudio de receptores de proteínas mediante el uso de redes complejas para el desarrollo y diseño de fármacos eficaces en patologías del sistema nervioso central. Memoria presentada por Manuel Quintín Escobar Cubiella Para optar al grado de Doctor en Farmacia Departamento de Química Orgánica Facultade de Farmacia Directores: Dr. Francisco Prado Prado Prof. Dr. Xerardo García Mera Dr. Humberto González Díaz Santiago de Compostela, Abril 2012 - 2 - - 3 - D. Xerardo García Mera, Prof. Titular y D. Francisco Javier Prado Prado, PDI Doctor Contratado por el Programa Ángeles Albariño, ambos del Departamento de Química Orgánica de la Universidad de Santiago de Compostela (USC), así como D. Humberto González Díaz, Doctor Contratado por el Departamento de Microbiología y Parasitología, Área de Parasitología, Facultad de Farmacia, USC. CERTIFICAN: Que la memoria titulada: “MODELOS BIOINFORMÁTICOS Y ESTUDIO DE RECEPTORES DE PROTEÍNAS MEDIANTE EL USO DE REDES COMPLEJAS PARA EL DESARROLLO Y DISEÑO DE FÁRMACOS EFICACES EN PATOLOGÍAS DEL SISTEMA NERVIOSO CENTRAL”, que para optar al grado de Doctor Farmacia presenta MANUEL QUINTÍN ESCOBAR CUBIELLA, ha sido realizada bajo nuestra dirección, en el Departamento de Química Orgánica de la Facultad de Farmacia de la Universidad de Santiago de Compostela. Y considerando que el trabajo constituye tema de Tesis Doctoral, autorizamos su presentación en la Universidad de Santiago de Compostela. Y para que conste, expedimos el presente certificado en Santiago de Compostela a nueve de abril del dos mil doce. ____________________________ _____________________________ Fdo.: Prof. Dr. Xerardo García Mera Fdo.: Dr. Francisco Javier Prado Prado ____________________________ Fdo.: Dr. Humberto González Díaz - 4 - - 5 - Castro soy y Castro he sido Asiento en firme Montaña Y a la Corona de España Con lealtad siempre he servido Armas, Escudo y Señal Castillo, Puente y Santa Ana Naos, Ballena y mar llana Son de Castro la Leal Divisa de la Ciudad de Castro-Urdiales. Siglo XIII-XIV. - 6 - - 7 - Agradecimientos Hace 20 años ya, que comenzó mi relación con la Universidad de Santiago de Compostela, y puedo afirmar, que desde el primer momento se afianzó en mi la ambición por el saber y el conocimiento. En primer lugar, tuve la enorme suerte de vivir mis experiencias universitarias en un entorno multicultural, el C.M.U. Gelmírez, donde me embadurné de amigos y compañeros geniales, de orígenes e inquietudes tan dispares con los que las tertulias y disputas dialécticas más alocadas se hacían sabrosamente interminables. En la Facultad de Farmacia fui descubriendo diversos aspectos de la profesión que aún sigo aplicando a día de hoy, pero especialmente me quedé enganchado de la asignatura de Química Farmacéutica, a mi entender, una de las disciplinas más completas de la carrera, y cuando ya había leído la Tesina y realizado los cursos de Doctorado, entro en contacto con la Industria Farmacéutica, teniendo el privilegio de haber aprendido y trabajado en grandes compañías, JanssenCilag, Merck Sharp & Dohme y los últimos 8 años en Novartis Pharmaceutica, y mira por donde mi camino se vuelve a cruzar con la USC a fin de rematar otro paso con la inquietud y ambición de no permitir que sea el último. Quiero dedicar esta memoria a las personas que con su esfuerzo y sacrificio han facilitado que me pudiera mover por el mundo adelante sin más preocupación que la de buscar mi sitio. Gracias a mis padres Manolo y Rebeca por su ayuda constante e incondicional y por el ejemplo que siempre he recibido de vosotros. A Cristina quiero agradecerle su comprensión por el tiempo robado y por el amor que me ha dedicado siempre, confiando ciegamente en mis posibilidades, además de por los dos hijos más guapos del mundo, Mikel y Helena, que algún día seguirán sus propias inquietudes. Espero saber estar a la altura…. A mi hermano Javier, primer Doctor de la Familia, porque me ha servido de acicate. A Xerardo García Mera, que desde el año 1994 me ha espoleado y prestado su modelo en la forma de entender la Ciencia. A Fran, por hacerme llegar a entender conceptos que ni sabía que existían, gracias por tu paciencia para conmigo. A Humber por su inestimable colaboración. - 8 - - 9 - Abreviaturas utilizadas LDA: Linear Discriminant Analysis, término que proviene del inglés: Análisis Discriminante Lineal. AChE: Acetilcolinesterasa. ANN: Artificial Neural Networks, término que proviene del inglés: Redes neuronales artificiales. 3D: Tridimensional. CM: Cadenas de Markov. DTP: Pares de fármaco-proteína con alta afinidad. nDTP: Pares de fármaco-proteína con nula afinidad. EA: Enfermedad de Alzheimer. EP: Enfermedad de Parkinson. FDA: Food and Drug Administration of USA; Administración de alimentos y medicamentos de los EE.UU. HTS: High-Throughput-Screening, término que proviene del inglés: evaluación de alta eficacia. IT: Índices topológicos LNN: Lineal Neural Network, término que proviene del inglés: Red neuronal lineal. MARCH-INSIDE (MI): Markov Chain Invariants for Network Simulation and Design PDB: Protein Data Bank; Banco de datos de proteínas. QSAR: Quantitative-Structure-Activity-Relationship, término que proviene del inglés: relación-cuantitativa-estructura-actividad. mt-QSAR: multi-target Quantitative-Structure-Activity-Relationship, término que proviene del inglés: relación-cuantitativa-estructura-actividad multi-target. QSPR: Quantitative-Structure-Property-Relationship, término que proviene del inglés: relación-cuantitativa-estructura-propiedad. QSTR: Quantitative-Structure-Toxicity-Relationship, término que proviene del inglés: relación-cuantitativa-estructura-toxicidad. SN: Sistema Nervioso. - 16 - La búsqueda y desarrollo de fármacos eficaces para el tratamiento de estas enfermedades ha generado grandes expectativas, debido a la relevancia que tienen sobre la economía de los sistemas sanitarios y la tremenda carga y desgaste que sufren familia y cuidadores. Por ello, la industria farmacéutica se ha volcado sobre estas patologías en las últimas tres décadas, pero las dificultades de realizar ensayos sobre el SN provoca que los gastos y tiempos de investigación se disparen, limitando de forma considerable la rentabilidad de los procesos tradicionales en el desarrollo de nuevos medicamentos8. Es en este apartado donde realiza sus aportaciones el diseño de fármacos, dedicando una parte del mismo al desarrollo de modelos matemáticos que permitan predecir propiedades de interés para una gran variedad de sistemas químicos incluyendo moléculas de bajo peso molecular, polímeros, biopolímeros, sistemas heterogéneos, formulaciones farmacéuticas, conglomerados de moléculas e iones, materiales, nanoestructuras y otros.9 Este tipo de predicciones tienen como objetivo fundamental complementar y evolucionar las técnicas de carácter experimental tradicionales, fundamentalmente colaborando en la obtención de nuevas moléculas activas con mayor probabilidad de éxito, con la ventaja que de ello se deriva en términos de ahorro de tiempo, recursos materiales y también en el refinamiento y reducción en el uso de animales de laboratorio.10-12 En la actualidad existen miles de compuestos químicos, ya sea de origen natural o de síntesis, y en su amplia mayoría aún no se les ha encontrado aplicaciones farmacológicas, agroquímicas, industriales o de algún otro tipo. Esto es consecuencia directa de la gran diferencia que existe entre la velocidad con que los nuevos compuestos son obtenidos y caracterizados en la mesa de laboratorio y la realización de los ensayos experimentales que permitan evaluar su potencial terapéutico. Cuando estos productos pretenden ser destinados al consumo humano en aras de conseguir mejoras terapéuticas, los procesos aún sufren una ralentización adicional por la cantidad de trámites administrativos y burocráticos destinados a garantizar que los fármacos lleguen al consumo humano con todas las garantías necesarias para los pacientes. Por otra parte existen determinadas patologías en las cuales los ensayos de laboratorio son de una gran complejidad en sí mismos, como ocurre en el caso de los compuestos con potencial actividad antiviral, o bien los ensayos son muy costosos tanto en términos de recursos materiales, humanos y de tiempo, como ocurre con los compuestos con potencial actividad neuroprotectora dirigidos al tratamiento de diferentes enfermedades degenerativas como EA y EP. - 17 - Este motivo ha obligado a la industria farmacéutica a cambiar sus estrategias de búsqueda enfocando sus esfuerzos hacia el desarrollo de métodos que racionalicen los sistemas de diseño y evaluación de nuevos fármacos ya en las primeras fases del descubrimiento de los mismos, produciéndose así una importante disminución del coste económico y temporal en el desarrollo de nuevos compuestos para su uso farmacéutico. Por ello, y en los últimos años, ha ido ganando una gran importancia el desarrollo de modelos capaces de predecir las rutas de descubrimiento con mayor probabilidad de éxito y que sirvan de guía al investigador en el desarrollo racional de fármacos con un importante ahorro de recursos. En dicho sentido, los estudios QSAR (Quantitative Structure-ActivityRelationships) son usados cada vez mas como herramientas para el descubrimiento molecular. Estos modelos QSAR pueden ser diseñados para que predigan la probabilidad de que un fármaco sea efectivo contra una enfermedad degenerativa determinada ya sea la enfermedad de Parkinson, Alzheimer o cualquier otra, actuando sobre una diana molecular específica. Para ello, primeramente deberán recopilarse de las bases de datos públicas los datos de actividad biológica de fármacos neuroprotectores con diana molecular conocida. Seguidamente, se calcularán determinados parámetros numéricos llamados Índices Topológicos (ITs) tanto de los fármacos como de sus dianas moleculares13 y posteriormente por análisis estadístico y/o Inteligencia Artificial (Redes Neuronales Artificiales) se buscarán los modelos QSAR. Dichos ITs describen únicamente la topología ó conectividad en fármacos y secuencias de sus dianas (proteínas y/o ADN/ARN) y, por ello, son versátiles y fáciles de usar. Además se pueden explorar grandes bases de datos con un importante ahorro de tiempo y recursos materiales14. Estos modelos duales QSAR Bioinformáticos + Quimioinformáticos son de interés para los grupos de investigación dedicados a la Química Farmacéutica. Se pueden usar estos modelos para predecir la probabilidad de actuar como dianas de fármacos a nuevas proteínas y/o ADN/ARN que se aíslen y que actúen dentro del desarrollo de las enfermedades demenciales, así como otras ya conocidas pero con función desconocida como diana de fármacos. En esta memoria presentamos de manera conjunta la revisión de modelos previos y trabajos específicos novedosos, en los que se han introducido nuevos índices numéricos utilizados para describir tanto la estructura molecular de fármacos como la estructura macromolecular de sus dianas o receptores (proteínas y/o ADN/ARN). Con - 18 - estos ITs hemos sido capaces de desarrollar nuevos modelos multiQSAR de gran interés por su doble función en la predicción de fármacos y sus dianas moleculares. Estos trabajos permitirán la introducción de nuevos conceptos teóricos y la evolución hacia modelos con posibles aplicaciones en la búsqueda de nuevos fármacos neuroprotectores útiles en el tratamiento de las enfermedades de Parkinson y Alzheimer y/o nuevas dianas moleculares para estos fármacos. Este tipo de investigación abarca un área general-básica en la que interactúan la Bioinformática y la Quimioinformática. 1.1 Metodología QSAR Esta metodología se basa en el uso de cálculos por ordenador y en las nuevas tecnologías de la informática las cuales pueden ser usadas tanto para pequeñas moléculas como para macromoléculas. Para moléculas pequeñas: 1. Estudios de relación cuantitativa estructura molecular-actividad farmacológica (QSAR) y de estructura molecular-propiedades toxicológicas y eco-toxicológicas incluyendo mutagenicidad y carcinogénesis (QSTR). 2. Predicción de propiedades químicas y fisicoquímicas de moléculas. Estudios de relación estructura molecular y propiedades de absorción, distribución, metabolismo y eliminación (ADME). 3. Predicción de mecanismos de acción biológica de moléculas y evaluación in silico de alta eficacia para grandes bases de datos (virtual HTS). Para macromoléculas: 4. Estudios de interacción fármaco-receptor (neuronas). 5. Bioinformática aplicada a estudios de relación secuencia-función y propiedades estructurales de ácidos nucleicos y proteínas. 6. Búsqueda de nuevas dianas terapéuticas y “sitio activo” a partir de datos de Genómica y/o Proteómica. 7. Búsqueda de biomarcadores para diagnóstico de enfermedades o como indicadores de contaminación. 8. Predicción de propiedades fisicoquímicas de polímeros sintéticos, biopolímeros, materiales y nano-estructuras. 9. Predicción, diseño y optimización de enzimas mutadas para procesos biotecnológicos. - 19 - 1.2. Desarrollo de la Metodología (Pasos a seguir) Entre las técnicas de regresión usadas, son de destacar las técnicas lineales debido a su sencillez. En ellas se intenta modelar la actividad biológica como una función lineal multivariada de los descriptores moleculares. Por otra parte, en etapas tempranas del descubrimiento molecular así como del estudio del mecanismo de acción de los fármacos, es suficiente tener una respuesta acerca de la probabilidad con que un fármaco tendrá la actividad o mecanismo bajo estudio, sin predecir el valor exacto. Particularmente, en este trabajo se utilizará el Análisis Discriminante Lineal (LDA) ya que posee la cualidad de ser simple permitiendo la clasificación de objetos en grupos predeterminados basándose en múltiples rasgos. En nuestro caso los objetos serán moléculas y los grupos el grado de actividad o un mecanismo de acción determinado. La estrategia general de trabajo en QSAR con LDA se puede dividir en una serie de pasos que son ilustrados gráficamente en la Figura 113, 14: 1. Recopilación de una serie de datos aleatoria, representativa y estratificada de moléculas con la actividad deseada y un grupo control que no posee la actividad bajo estudio. 2. Selección de los descriptores moleculares a utilizar. 3. Cálculo, mediante un programa computacional, de los descriptores moleculares seleccionados. 4. Utilización de los descriptores calculados a las moléculas recopiladas (serie de entrenamiento) para determinar modelos QSAR en un programa de cálculo estadístico. 5. Validación de los modelos QSAR contrastando la actividad predicha a las moléculas recopiladas (serie de predicción) con su actividad experimental. 6. Uso de los modelos encontrados para predecir la actividad a moléculas no ensayadas con anterioridad. - 20 - Figura1: Esquema general de trabajo con técnicas QSAR. - 21 - 2. Los Descriptores Moleculares El número de moléculas que puede ser obtenido por síntesis orgánica es tan elevado que la probabilidad de seleccionar al azar una molécula que presente la actividad biológica deseada es prácticamente nula. Como consecuencia, ha surgido un gran interés en los métodos teóricos que relacionan la estructura molecular con la actividad biológica, entre ellos especialmente el QSAR. Uno de los pilares para el desarrollo del QSAR lo constituye el proceso de codificar la estructura química mediante descriptores moleculares.15 Más estrictamente, un descriptor molecular es el resultado final de una lógica y de un procedimiento matemático que transforma la información química codificada dentro de una representación simbólica de una molécula en un número útil o el resultado de algún experimento estandarizado. Estas variables pueden ser teóricas o experimentales, pueden describir a la molécula como un todo (descriptores globales) o solo representar un fragmento presente en ella (descriptores fragmentos). Se han definido en la actualidad más de 4000 descriptores moleculares diferentes, véase por ejemplo el programa de cálculo DRAGON y el HANDBOOK de descriptores moleculares recopilados por sus autores (Figura 2).16 En la actualidad y en relación con la constante definición de nuevos índices estructurales o descriptores moleculares no se detecta la “explosión” vista en pasados años, observándose cierta tendencia a una aplicación más diversificada e intensiva de los mismos. No obstante, el enorme número de propiedades a estudiar condiciona la continua introducción de nuevos índices basados en otros métodos, con la intención de que los químico-farmacéuticos posean un “arsenal” de descriptores moleculares lo más completo posible. Este contexto ha propiciado que algunos investigadores se hayan dado a la tarea de crear fórmulas matemáticas de los descriptores moleculares que ofrezcan un cuadro unificado de los mismos, para facilitar su sistematización y estudio.17-19. Estas fórmulas pueden, además, indicar la dirección de búsqueda para nuevos descriptores moleculares. En este trabajo se han usado dos programas bien conocidos para el cálculo de descriptores moleculares: el software MARCH INSIDE (MI) y el DRAGON. - 22 - Figura 2. Interface del programa DRAGON que calcula más de 4000 descriptores moleculares agrupados en 18 familias diferentes. 2.1. Descriptores del DRAGON El programa DRAGON ofrece la posibilidad de calcular un gran número de descriptores moleculares agrupados en diferentes familias. A su vez, la lista de descriptores proporcionados puede ser organizada como cerodimensionales (0D), unidimensionales (1D), bidimensionales (2D) y tridimensionales (3D). En nuestro caso y con el objetivo de simplificar la descripción utilizaremos esta última clasificación.20 Los descriptores calculados en este trabajo son obtenidos con la aplicación de este programa y son cantidades teórico-definidas no habiendo utilizado en ningún caso descriptores experimentales. Descriptores 0D: describen solamente la constitución de la molécula, pero no aportan información concerniente a la conformación ni al tipo de conectividad presente. Los más simples son, entre otros, el número de átomos de un determinado tipo, el número de enlaces y el peso molecular. Descriptores 1D: describen fragmentos de las moléculas formados por el agrupamiento de sus átomos constituyentes. Descriptores 2D: utilizan una función de auto-correlación bidimensional que contiene la topología del grafo y, además, representa la distribución de una propiedad atómica - 23 - determinada en la molécula. La propiedad atómica con la que se pesa/pondera el descriptor considera los átomos presentes en la molécula a través de la electronegatividad, la masa atómica, la polarizabilidad atómica, el estado electrotopológico o el volumen de Van der Waals, con lo cual se pueden seleccionar aquellos átomos que proporcionen mayor peso a la variable considerara. Estos descriptores tienen en cuenta las interacciones inter/intramoleculares. Descriptores 3D: esta clase tiene en cuenta los aspectos conformacionales de la estructura molecular, considerando de esta manera las propiedades estereoquímicas de las moléculas. Para su cálculo se utilizan estructuras moleculares previamente optimizadas con métodos convenientes, tales como el Método de campos de fuerza de la mecánica molecular MM+, en combinación con métodos derivados de la Mecánica cuántica, sean ab initio o Métodos de la teoría semiempírica de orbitales moleculares. Entre estos descriptores citamos las cargas atómicas, la energía del orbital molecular más alto ocupado y la energía del orbital molecular más bajo desocupado, entre otros. Un descriptor debe cumplir con un conjunto de características tales como: i. Fácil cálculo. ii. Invarianza respecto de la traslación y la rotación. iii. Invarianza respecto a la numeración de los átomos. iv. Buena correlación con la propiedad estudiada. v. Bajo grado de correlación con otros descriptores. 2.2. Descripción del método MARCH-INSIDE. Cadenas de Markov (CM) es el nombre de una teoría o tipo de modelo matemático definido por Markov.21 En nuestro trabajo hemos utilizado especialmente el método MARCH-INSIDE (del inglés: Markov Chain Invariants for Network Simulation and Design), el cual emplea las CM para calcular descriptores moleculares mediante una aproximación sencilla a fenómenos tales como 22-27: i. distribución de electrones de valencia alrededor de los átomos de una molécula. ii. propagación de una vibración en una cadena de RNA. iii. propagación de interacciones electrostáticas superficiales en una proteína viral o en la estructura plegada 3D de una enzima. iv. paso átomo por átomo de un fármaco desde el plasma a un tejido. v.interacción paso por paso de un fármaco con su receptor. - 24 - - 25 - 3. Estudios Teóricos Fármaco – Proteína La predicción rápida y precisa de las interacciones entre los fármacos y proteínas es una pieza clave en la combinación de la bioinformática y la investigación del proteoma hacia el descubrimiento de fármacos. Por lo tanto, hay un fuerte incentivo para desarrollar nuevos métodos capaces de detectar estas posibles interacciones fármaco-proteína de manera eficiente.28 En este sentido, los gráficos y la teoría de redes complejas pueden jugar un papel importante en las diferentes etapas del proceso de modelado con diferentes grados de organización de la materia.29-36 3.1. Etapas de los estudios Interacción Fármaco-Proteína. En una primera etapa, podemos utilizar los gráficos moleculares no sólo para representar y calcular los parámetros estructurales de los fármacos, los ITs, sino también para estimar los parámetros físico-químicos sobre la base de un método gráfico.37 En un nivel superior, podemos utilizar los gráficos para representar la estructura de los fármacos-proteínas y calcular los ITs característicos y/o los parámetros físicoquímicos de la estructura de las proteínas o las redes de interacciones entre proteínas, véase por ejemplo las obras de Giuliani,38-44 o las revisiones publicadas en los últimos años.45 A continuación, se pueden utilizar los ITs y/o los parámetros físico-químicos como entradas para la búsqueda de clasificadores lineales o no lineales capaces de predecir la red -como las estructuras moleculares que presentan o no una propiedad de interés, ver por ejemplo las obras de Caballero y Fernández et al.46-49con aplicaciones tanto al campo de los medicamentos y las proteínas o las obras de Zbilut et al.50,51 En particular, utilizando los parámetros del fármaco y de la proteína se puede discriminar entre pares fármaco-proteína con alta afinidad (DTPs) y pares fármacoproteína con nula afinidad (nDTPs). El método QSAR puede convertirse en una herramienta muy útil en este contexto para reducir sustancialmente el tiempo y los recursos que consumen los experimentos. En una última etapa, la predicción de todos los posibles DTPs/nDTPs de la base de datos forman la red compleja de fármacos y proteínas. Por ejemplo, Yildirim y Goh et al.52 han construido un grafo bipartito compuesto por la base de datos de la FDA, en - 32 - - 33 - 4.3. Revisión de estudios teóricos y bioinformáticos de inhibidores de la acetilcolinesterasa. Manuel Escobar, Franco Fernández, Xerardo García Mera y Francisco Javier Prado Prado Current Bioinformatics, 2012, in press. La enfermedad de Alzheimer (EA) es la enfermedad del siglo XXI, y no sólo es compleja y difícil de tratar sino que hasta ahora es incurable. Tampoco sabemos con certeza que podemos hacer para prevenirla. Es por ello que los tratamientos actuales se centran en múltiples aspectos incluyendo la ayuda a las personas a mantener la función mental, los síntomas conductuales y frenar, retrasar o prevenir la enfermedad. En la actualidad existen cuatro medicamentos aprobados por los EE.UU, Food and Drug Administration (FDA) para el tratamiento de la EA, donepezilo, rivastigmina y galantamina se usan para el tratamiento del Alzheimer en fases de leve a moderada y la memantina se utiliza para tratar los niveles moderado a severo. Estos medicamentos actúan regulando los neurotransmisores (las sustancias químicas que transmiten mensajes entre las neuronas). El tratamiento de la EA por los precursores de acetilcolina y los agonistas colinérgicos se ha mostrado ineficaz e incluso ha originados graves efectos secundarios. La hidrólisis de la ACh por la AChE causa la terminación de la neurotransmisión colinérgica. Por lo tanto, los compuestos que inhiban significativamente la AChE podrían aumentar los niveles de ACh y paliar la EA. Sin embargo, estos medicamentos no cambian el proceso de la enfermedad y pueden ayudar unos pocos meses o incluso unos pocos años. En este sentido, la metodología QSAR podría desempeñar un papel importante en el estudio de estos inhibidores de la acetilcolinesterasa (ACE). Los modelos QSAR pueden contribuir a guiar la síntesis de la AChE. En este trabajo hemos revisado diferentes estudios de bioinformática y estudios teóricos de inhibidores de la AChE, estudios de diseño y cálculo de una serie muy grande y heterogénea de los inhibidores de este enzima. En primer lugar, revisamos los métodos 2D QSAR, 3D QSAR, CoMFA, CoMSIA y de docking con diferentes compuestos para averiguar los requisitos estructurales. A continuación, revisamos los estudios QSAR utilizando el método de LDA con el fin de comprender las exigencias estructurales esenciales para la unión con el receptor de dichos inhibidores de la AChE, utilizando descriptores de ModesLab de una base de datos de 10.000 fármacos diferentes del servidor ChemBL. - 34 - - 35 - 5. Objetivos 5.1. Objetivos Generales: I. Desarrollar nuevos modelos QSAR aplicables a la predicción de la actividad biológica de compuestos contra una única diana o múltiples dianas, modelos QSAR multi-target (mt-QSAR), de interés en química farmacéutica, microbiología, y parasitología. II. Desarrollar nuevas metodologías haciendo uso de varios programas de cálculo de descriptores moleculares para la construcción de Redes Complejas de compuestos útiles en estudios de Bioinformática a partir de modelos QSAR ó mt-QSAR. 5.2. Objetivos específicos: I. Desarrollar modelos QSAR para la predicción de inhibidores de la β-secretasa. II. Desarrollar modelos QSAR para la predicción de inhibidores de la acetilcolinestarasa. III. Desarrollar modelos mt-QSAR para la predicción de inhibidores de la acetilcolinestarasa. IV. Desarrollar nuevas metodologías para generar nuevos modelos que permitan estudiar las interacciones fármaco-proteína. V. Desarrollar una metodología QSAR/QSPR usando nuevos índices topológicos para construir redes de interacción fármaco-proteína. VI. Desarrollar un nuevo método para cuantificar numéricamente la calidad de las conexiones entre vértices basándose en los índices Markov para redes fármacoproteína. - 36 - - 37 - En este punto se presentarán todos los resultados obtenidos en forma de artículos de investigación ya publicados por el autor. Los 4 artículos presentados están agrupados de acuerdo al objetivo específico que cumplimentan. Para cada artículo se presenta una breve sección explicativa en Español de su importancia y los resultados alcanzados. En el apartado “9. Anexos (Publicaciones)” de esta Tesis se adjuntan las publicaciones correspondientes en el idioma en que fueron publicadas. 6. Resultados y Discusión - 38 - - 39 - 6.1. 2D MI-DRAGON: un nuevo modelo para estudiar las interacciones proteínaligando y el estudio teórico-experimental de la red fármaco-proteína obtenida de la base de datos de la FDA de EE.UU., estudio de oxoisoaporfinas como inhibidores de la MAO-A y proteínas de parásitos en humanos. Francisco Prado-Prado, Xerardo García-Mera, Manuel Escobar, Eduardo SobarzoSánchez, Matilde Yañez, Pablo Riera-Fernández and Humberto González-Díaz. European Journal of Medicinal Chemistry, 2011, 46 (12), 5838-5851. Hay muchos posibles pares de fármacos-proteínas que pueden tener lugar o no (DTP/nDTPs) entre las drogas con alta afinidad/no afinidad con proteínas diferentes. Por este motivo resulta costoso en términos de tiempo y recursos, por ejemplo, la determinación de todas las posibles interacciones entre ligandos de proteínas para un solo medicamento. En este aspecto, podemos utilizar el QSAR para llevar a cabo la predicción de DTP. Desafortunadamente, casi todos los modelos QSAR predicen la actividad contra un solo objetivo, proteína o diana. Para solucionar este problema se puede desarrollar un multi-target QSAR (mt-QSAR). En este trabajo, presentamos la técnica 2D MI-DRAGON un nuevo predictor de DTP basado en dos programas de software diferentes conocidos. Utilizamos el software MI para el cálculo de parámetros estructurales en 3D para las proteínas y el software DRAGON, uno de los más completos y basado en cálculos con más de 1600 parámetros, para calcular dichos descriptores moleculares en 2D mostrando todos los fármacos que tienen interacción conocida con proteínas presentes en la base de datos de la FDA. Un modo de desarrollar este tipo de ensayos en mt-QSAR consiste en incorporar en las ecuaciones QSAR parámetros de la estructura de las dianas (proteínas, DNA, RNA, etc..) añadidos a los parámetros estructurales de los fármacos presentes en el clásico QSAR. Ambas clases de parámetros fueron utilizados como inputs en diferentes algoritmos de redes neuronales artificiales para buscar un modelo no lineal mt-QSAR muy preciso. El mejor modelo ANN encontrado es un perceptrón multicapa (MLP) cuyo perfil es MLP 21:2131-1:1, el cual clasifica correctamente 303 de 339 DTP (sensibilidad = 89,38%) y 480 de 510 nDTPs (especificidad = 94,12%), que corresponde al promedio de la serie de entrenamiento = 92,23%. La validación del modelo se llevó a cabo por medio de una serie externa de predicción con una sensibilidad = 92,18% (625/678 DTP) y con una especificidad = 90,12% (730/780 nDTPs) y un promedio = 91,06%. - 40 - El modelo 2D MI-DRAGON ofrece una buena oportunidad para calcular de forma rápida todos los DTP posibles de un fármaco lo que nos permite reconstruir grandes redes complejas de DTP. Por ejemplo, hemos reconstruido con la base de datos de la FDA una red compleja con 855 nodos de 519 medicamentos + 336 objetivos. Hemos predicho una red compleja con topología similar (valores observados y previstos de media distancia es igual a 6,7 frente a 6,6. Esta red compleja puede ser utilizada para explorar grandes bases de datos de DTP con el fin de descubrir los nuevos medicamentos y / o proteínas. Finalmente, se ilustra con un estudio teórico-experimental el uso práctico del modelo 2D MI-DRAGON. En este trabajo se realiza la predicción, síntesis y evaluación farmacológica de 10 oxoisoaporfiinas diferentes con actividad inhibitoria de la MAO-A. El compuesto más activo OXO5 presentó una IC50 = 0,83 μM, notablemente mejor que la clorgilina que se tomó como fármaco control. Es posible buscar excelentes predictores de DTP utilizando como entrada parámetros estructurales de los fármacos y proteínas calculados con diferentes programas y en combinación con modelos ANN. El modelo 2D MI-DRAGON, basado en los parámetros estructurales de los fármacos calculada con el DRAGON y los parámetros de proteínas calculado con el MI, predice correctamente los DTPs de 500 fármacos diferentes aprobados por la FDA con una precisión mayor del 90%. El modelo 2D MI-DRAGON también es útil para construir y desarrollar nuevas redes complejas computacionalmente construidas, las cuales ofrecen una alternativa para descubrir nuevos medicamentos y proteínas y explorar la selectividad y la toxicidad de los fármacos. En este trabajo estas conclusiones se ejemplifican a través del estudio experimental y teórico de las nuevas isoaporfinas que presentan actividad inhibitoria de la MAO-A. - 41 - 6.2. 3D MI-DRAGON: nuevo modelo para la reconstrucción de la red fármacoproteína de la base de datos de la FDA y estudios teóricos experimentales de los derivados de la rasagilina como inhibidores de la acetilcolinesterasa. Francisco Prado-Prado, Xerardo García-Mera, Manuel Escobar, Nerea Alonso, Olga Caamaño,1 Matilde Yañez and Humberto González-Díaz. European Journal of Organic Chemistry, 2012. Submitted Las enfermedades neurodegenerativas se han incrementado de manera notable en los últimos años. Muchos de los fármacos que se utilizan en el tratamiento de dichas enfermedades presentan características específicas estructurales en 3D. Una proteína importante en este sentido es la acetilcolinesterasa, que es la diana de muchos fármacos usados en la Enfermedad de Alzheimer. En consecuencia, la predicción de las interacciones fármaco-proteína (DTP/nDTPs) entre nuevos fármacos candidatos con estructura 3D específica y proteínas adquiere gran relevancia. Por ello, podemos utilizar técnicas mt-QSAR para realizar la predicción de DTPs. Desafortunadamente, muchos de los modelos QSAR previamente desarrollados para predecir DTPs toman en consideración únicamente la información estructural 2D y codifican la actividad frente a una sola diana o proteína. Para contribuir a la resolución de este problema, hemos desarrollado un modelo 3D multi-target QSAR (3D mt-QSAR). En este trabajo, se introduce la técnica de 3D MI-DRAGON un nuevo modelo que predice los DTPs basándose en la utilización de dos clases de software conocidos. Así utilizamos el software MI y DRAGON 3D para el cálculo de los parámetros estructurales de los fármacos y las proteínas, respectivamente. Ambas clases de parámetros 3D se utilizan como materia prima para entrenar la ANN, utilizando los datos como referencia para construir la red compleja formada por todos los DTPs encontrados como fármacosproteínas que han sido aprobadas por la FDA. El conjunto de datos se ha descargado del Drug Bank. El mejor modelo 3D mt-QSAR encontrado es un ANN de tipo Multi-Layer Perceptron (MLP) con el perfil MLP 37:37-24-1:1, el cual clasifica correctamente 274 de 321 DTP (sensibilidad = 85,35%) y 1041 de 1190 nDTPs (especificidad = 87,48%), que corresponde al promedio = 87,03%. La validación del modelo se ha llevado a cabo con una serie de predicción externa con sensibilidad = 84,16% (542/644 DTP). Especificidad = 87,51% (2039/2330 nDTPs) y Promedio = 86,78%. Las nuevas redes complejas de los DTPs han sido reconstruidas a partir de la base de datos de la FDA y - 48 - - 49 - 1. Olson, R.E. Annu. Rep. Med. Chem. 2000, 35, 31-40. 2. Woodgett, J.R. EMBO J. 1990, 9 (8), 2431-2438. 3. Troussard, A.A; Tan, C., Yoganathan, T.N. Dedhar S. Molecular Cell Biology. 1999, 19, 7420-7427. 4. Turenne, G.A.; Price, B.D. BMC Cell. Biol. 2001, 2, 12-21. 5. Hoeflich, K.P.; Luo, J.; Rubie, E.A.; Tasao, M.S.; Jin, O.; Woodgett, J.R. Nature. 2000, 406, 86-90. 6. Hooper, C.; Killick, R.; Lovestone, S. J. Neurochem. 2008, 104, 1433-1439. 7. MacAulay, K.; Doble, B.W.; Patel, S.; Hansotia, T.; Sinclair, E.M.; Drucker, D.J. et al. Cell Metab. 2007, 6, 329-337. 8. Droucheau, E.; Primot, A.; Thomas, V.; Mattei, D.; Knockaert, M.; Richardson, C.; et al. Biochim Biophys. Acta. 2004, 1697, 181-196. 9. Kubinyi, H.; Taylor, J.; Ramdsen, C. Quantitative Drug Design, in Comprehensive Medicinal Chemistry. Ed. C. Hansch. Pergamon. 1990, 4, 589643. 10. Lutz, M.W.; Menius, J. A.; Laskody, R.G.; Domanico, P.L.; Goetz, A.G.; Saussy, D. L.; Rimele, T. Network Science. 1996, 2, Issue 9, September. 11. Loew, G.H.; Villar, H.O.; Alkorta, Y. Pharm. Res. 1993, 10, 475-486. 12. Wess, G. Drug Discovery Today. 1996, 1, 529-535. 13. Kier, L.B.; Hall, L.H. Topological Indices and Related Descriptors in QSAR and QSPR. Gordon and Breach, Amsterdam, 1999. 14. Todeschini, R.; Consonni, V. Handbook of Molecular Descriptors. Mannhold, R.; Kubinyi, H.; Timmermann, H. Ed., Wiley-VCH: Weinheim, 2000. 15. Devillers, J.; Balaban, A.T. Topological Indices and Related Descriptors in QSAR and Drug Design, Amsterdan, 2000, pp 3-41. 8. Referencias Bibliográficas - 50 - 16. Todeschini, R.; Consonni, V. Handbook of Molecular Descriptors. Wiley VCH, Weinheim, Germany. 2000. 17. Estrada, E.; Chem. Phys. Lett .2001, 336, 9890-9895. 18. Kier, L.B.; Hall, L.H. Topological Indices and Related Descriptors in QSAR and QSPR. Gordon and Breach Sci. Pub.: Amsterdam, 1999, pp 455-489. 19. Wiener, H. J. Am. Chem. Soc. 1947, 69, 17-20. 20. Helguera, A.M.; Combes, R.D.; González, M.P.; Cordeiro, M.N.D.S.. Curr. Top. Med. Chem. 2008, 8 (18), 1628-1655. 21. Markov, A.A. Bull. Soc. Phys. Math. Kasan 1906, 15, 155-165. 22. Bharucha-Reid, A.T. Elements of Theory of Markov Process on the Application, McGraw-Hill Series in Probability and Statistic, McGraw-Hill Book Company, New York. 1960, 167-434. 23. Freund, J.A.; Poschel, T. Eds. Stochastic Processes in Physics, Chemistry, and Biology. In: Lect. Notes Phys. Springer-Verlag, Berlin, Germany 2000. 24. González-Díaz, H.; Prado-Prado, F.; Ubeira, F.M. Curr. Top. Med. Chem. 2008, 8 (18), 1676-90. 25. González-Díaz, H.; Duardo-Sanchez, A.; Ubeira, F.M.; Prado-Prado, F.; PérezMontoto, L.G.; Concu, R.; Podda, G.; Shen, B. Curr. Drug. Metab. 2010 May; 11 (4), 379-406. 26. González-Díaz, H.; González-Díaz, Y.; Santana, L.; Ubeira, F.M.; Uriarte, E. Proteomics. 2008, Feb, 8 (4), 750-778. 27. González-Díaz, H.; Vilar, S; Santana, L.; Uriarte, E. Curr. Top. Med. Chem. 2007, 7 (10), 1015-1029. 28. Yamanishi, Y.; Araki, M.; Gutteridge, A.; Honda, W.; Kanehisa, M., Bioinformatics 2008, 24 (13), i232-240. 29. Giuliani, A. BMC Genomics. 2010, 11 Suppl 1, S2. 30. Dhar, P. K.; Giuliani, A. Sys. Synth. Biol. 2010, 4, (1), 7-13. 31. Bornholdt, S.; Schuster, H. G. Handbook of Graphs and Complex Networks: From the Genome to the Internet. WILEY-VCH GmbH & CO. KGa. Wheinheim, 2003. 32. Estrada, E. Proteomics 2006, 6 (1), 35-40. 33. Estrada, E. J. Proteome Res. 2006, 5 (9), 2177-2184. 34. Réka, A.; Barabasi, A.L. Rev. Mod. Phys. 2002, 74 (1), 47-97. 35. Barabasi, A. L.; Oltvai, Z. N. Nat. Rev. Genet. 2004, 5 (2), 101-113. - 51 - 36. Barabasi, A. L. N. Engl. J. Med. 2007, 357 (4), 404-407. 37. González-Díaz, H.; Vilar, S.; Santana, L.; Uriarte, E. Curr. Top. Med. Chem. 2007, 7 (10), 1025-1039. 38. Giuliani, A.; Di Paola, L.; Setola, R. Curr. Proteomics 2009, 6 (4), 235-245. 39. Krishnan, A.; Zbilut, J. P.; Tomita, M.; Giuliani, A. Curr. Protein Pept. Sci. 2008, 9 (1), 28-38. 40. Krishnan, A.; Giuliani, A.; Zbilut, J. P.; Tomita, M. PLoS ONE. 2008, 3 (5), e2149. 41. Palumbo, M. C.; Colosimo, A.; Giuliani, A.; Farina, L. FEBS Lett. 2007, 581 (13), 2485-2489. 42. Krishnan, A.; Giuliani, A.; Zbilut, J. P.; Tomita, M. J. Proteome Res. 2007, 6 (10), 3924-3934. 43. Krishnan, A.; Giuliani, A.; Tomita, M. PLoS ONE. 2007, 2 (6), e562. 44. Tun, K.; Dhar, P. K.; Palumbo, M. C.; Giuliani, A. BMC Bioinformatics. 2006, 7-24. 45. Vilar, S.; González-Díaz, H.; Santana, L.; Uriarte, E. J. Theor. Biol. 2009, 261 (3), 449-458. 46. Caballero, J.; Fernández, M. Curr. Top. Med. Chem. 2008, 8 (18), 1580-1605. 47. Fernández, L.; Caballero, J.; Abreu, J. I.; Fernández, M. Proteins 2007, 67, 834– 852. 48. Fernández, M.; Caballero, J.; Fernández, L.; Abreu, J. I. J. Mol. Graph. Model. 2007, 26 (4), 748-759. 49. Fernández, M.; Caballero, F.; Fernández, L.; Abreu, J. I.; Acosta, G. Proteins 2008, 70 (1), 167-175. 50. Zbilut, J. P.; Giuliani, A.; Colosimo, A.; Mitchell, J. C.; Colafranceschi, M.; Marwan, N.; Webber, C. L., Jr.; Uversky, V. N. J. Proteome Res. 2004, 3 (6), 1243-1253. 51. Zbilut, J. P.; Colosimo, A.; Conti, F.; Colafranceschi, M.; Manetti, C.; Valerio, M.; Webber, C. L., Jr.; Giuliani, A. Biophys. J. 2003, 85 (6), 3544-3557. 52. Yildirim, M. A.; Goh, K. I.; Cusick, M. E.; Barabasi, A. L.; Vidal, M. Nat. Biotechnol. 2007, 25 (10), 1119-1126. 53. Kirchmair, J.; Markt, P.; Distinto, S.; Schuster, D.; Spitzer, G.M.; Liedl, K.R.; Langer, T.; Wolber, G. J. Med. Chem. 2008, 51, 7021-7040. - 52 - 54. Knox, C.; Law, V.; Jewison, T.; Liu, P.; Ly, S.; Frolkis, A.; Pon, A.; Banco, K.; Mak, C.; Neveu, V.; Djoumbou, Y.; Eisner, R.; Guo, A.C.; Wishart, D.S. Nucleic Acids Res. 2011, 39, 1035-1041. 55. Wishart, D.S.; Knox, C.; Guo, A.C.; Cheng, D.; Shrivastava, S.; Tzur, D.; Gautam, B.; Hassanali, M. Nucleic Acids Res. 2008, 36, 901-906. 56. Wishart, D.S.; Knox, C.; Guo, A.C.; Shrivastava, S.; Hassanali, M.; Stothard, P.; Chang, Z.; Woolsey, J. Nucleic Acids Res. 2006, 34, 668-672. - 53 - A continuación se presentan las diferentes publicaciones que se recogen en la Tesis siguiendo el orden de presentación de las mismas establecido en la memoria. 9. Anexos (Publicaciones) Current Bioinformatics, 2011, 6, 3-15 3 1574-8936/11 $58.00+.00 © 2011 Bentham Science Publishers Ltd Review of Bioinformatics and QSAR Studies of -Secretase Inhibitors Francisco Prado-Prado*, Manuel Escobar-Cubiella and Xerardo García-Mera Department of Organic Chemistry, University of Santiago de Compostela, Spain Abstract: Alzheimer disease (ADa) is the most common form of senile dementia, and it is characterized pathologically by decreased brain mass. An important problem to inhibiting -secretase, is to cross the blood-brain barrier (BBB) using drugs not derived from proteins and thus more efficient to design drugs to treat Alzheimer's disease. In this sense, quantitative structure-activity relationships (QSAR) could play an important role in studying these -secretase inhibitors. QSAR models are necessary in order to guide the -secretase synthesis. In the present work, we firstly revised two servers like ChEMBL or PDB to obtain databases of -secretase inhibitors. Next, we review previous works based on 2D-QSAR, 3DQSAR, CoMFA, CoMSIA and Docking techniques, which studied different compounds to find out the structural requirements. Last, we carried out new QSAR studies using Artificial Neural Network (ANN) method and the software ModesLab in order to understand the essential structural requirement for binding with receptor for -secretase inhibitors. Keywords: QSAR, CoMSIA, COMFA, docking, topological indices, -secretase inhibitors, alzheimer's disease (AD). 1. INTRODUCTION Alzheimer disease (ADa) is the most common form of senile dementia, and it is characterized pathologically by decreased brain mass, extracellular senile plaques, and intracellular neurofibrillary tingles [1]. The major risk factor for AD is age, and it affects more than a third of all people that reach 85 years of age. Although the pathogenesis of AD is still controversial, rare human mutations that lead to familial early onset AD are found in genes that lead to increased expression or processing of the amyloid precursor protein (APP) into amyloid  peptide (A), which is the major component of senile plaques [2-4]. The so-called “amyloid hypothesis” for AD pathogenesis is supported by transgenic mouse models that overexpress mutant APP and genes involved in APP processing, which lead to the production of senile plaques and cognitive impairment [5-8] (see Fig. 1). Prior to its identification, numerous studies were undertaken to define the characteristics of -secretase activity. Although the majority of body tissues exhibit -secretase activity [9], highest activity levels were observed in neural tissue and neuronal cell lines [10]. -Secretase also called BACE1 (-site of APP cleaving enzyme) (see Fig. 2), is an aspartic-acid protease important in the pathogenesis - of Alzheimer's disease, and in the formation of myelin sheaths in peripheral nerve cells [5-8]. The main problem inhibiting -secretase, is to cross the blood-brain barrier (BBB) using drugs not derived from proteins and thus more efficient to design drugs to treat Alzheimer's disease. There are two barriers in the brain: the blood-cerebrospinal fluid barrier and the BBB. The BBB consists of the endothelial cells that form the brain capillaries. It is restrictive to penetration of molecules, owing to tight junctions, lack of *Address correspondence to this autor at the Department of Organic Chemistry, University of Santiago de Compostela, Spain; Tel: + 881814940; Fax: 34-981-594912; E-mail: [email protected]. fenestrations, negative surface polarity, and the high level of efflux transporters [11]. In this part, Chemoinformatics and Bioinformatics methods may play an important role in the study of BACE1 inhibitors, Quantitative Structure-Activity Relationships (QSAR) studies are used as predictive tools for the molecular development [12, 13]. Up to today, there are near 1600 molecular descriptors that, in principle, can be generalized and used to solve the former problem [14]. Many of these indices are known as molecular Topological Indices (TIs) or simply invariants of a molecular graph. Unfortunately, QSAR studies are generally based on databases considering only structurally parent compounds acting against one single microbial species. In a recent review, our group have discussed recent advances in the field [15]. In addition to QSAR, Bioinformatics and Chemoinformatics methods useful to study -secretase may include techniques like Comparative Molecular Field Analysis (CoMFA), drug-target Docking, Sequence Alignment (SA) or other methods. In a recent, preliminary review in the field published in Proteomics in 2008 was discussed the use of these methods but only from the point of view of proteins [16]. Almost all QSAR techniques are based on the use of molecular descriptors, which are numerical series that codify useful chemical information and enable correlations between statistical and biological properties [17, 18]. On the other hand, QSAR models can be used to explore the relationships between the structural spaces of compounds as inhibitors for specific enzymes, such as MAO inhibitors [19], HIV-1 integrase inhibitors [20], and/or protease inhibitors [21] or tyrosinase inhibitors [22-24]. In fact, recently, the field has moved from small molecules to proteins and other systems. See for instance, recent review issues in Current Topics in Medicinal Chemistry [25-34], Current Proteomics [35-41], Current Drug Metabolism [42-50] and Current Pharmaceutical Design [51-60]. In the present work, we firstly revised servers like ChEMBL (http://www.ebi.ac.uk/ChEMBLdb/) or PDB (http://www.- 4 Current Bioinformatics, 2011, Vol. 6, No. 1 Prado-Prado et al. pdb.org/pdb/home/home.do) to obtain databases of - secretase inhibitors. Next, we review previous works based on 2D-QSAR, 3D-QSAR, CoMFA, CoMSIA and Docking techniques, which studied different compounds to find out the structural requirements. Last, we carried out new QSAR studies using Artificial Neural Network (ANN) method and the software ModesLab [61] in order to understand the essential structural requirement for binding with receptor for - secretase inhibitors. 2. THEORETICAL STUDIES FOR -SECRETASE INHIBITORS In this section we updated the contents presented in our recent review published in Current Drugs Metabolim [49]. The high number of possible candidates to -secretase inhibitors creates the necessity of Quantitative StructureActivity Relationship models in order to guide the - secretase inhibitor synthesis. In this work, we revised diffeFig. (1). APP metabolism by the secretase enzymes. Fig. (2). 3D of -secretase. Review of Bioinformatics and QSAR Studies Current Bioinformatics, 2011, Vol. 6, No. 1 5 rent computational studies for a very large and heterogeneous series of -secretase. First, we revised databases for bioinformatics studies of secretase for QSAR studies with conceptual parameters. Next, using method of regression analysis; and QSAR studies in order to understand the essential structural requirement for binding with receptor. Next, we review 3D QSAR, CoMFA, CoMSIA and Docking with different compound to find out the structural requirements for -secretase inhibitors. 2.1. Databases for Bioinformatics Studies of Secretase 2.1.1. ChEMBL Dataserver of -Seceretase Inhibitors for QSAR Studies ChEMBL [62] is a database of bioactive drug-like small molecules, it contains 2-D structures, calculated properties (e.g. logP, Molecular Weight, Lipinski Parameters, etc.) and abstracted bioactivities (e.g. binding constants, pharmacology and ADMET data) (see Fig. 3). In this review, the search for -secretase inhibitors by ChEMBL is important for the development of new theoretical studies and QSAR models or docking studies. The -secretase database attempt to normalise the bioactivities into a uniform set of end-points and units where possible, and also to tag the links between a molecular target and a published assay with a set of varying confidence levels. The data is abstracted and curated from the primary scientific literature, and cover a significant fraction of the SAR and discovery of modern drugs. Additional data on clinical progress of compounds is being integrated into ChEMBL at the current time. ChEMBL search by: • Search target data via keyword, protein sequence search (BLAST), or by navigating the target classification hierarchy. • Search compound data with lists of keywords, SMILES strings, or compound identifiers. Substructure and similarity searching functionality is available also. • Search assay data via keyword search using the main search bar. The ChEMBL database (ChEMBLdb) contains medicinal chemistry bioassay data, integrated from a wide variety of sources (the literature, deposited data sets, other bioassay databases). Subsets of ChEMBLdb, relating to particular target classes, or disease areas, are exported to smaller databases like -seceretase inhibitors. ChEMBL target search results for secretase enzymes was 11 hits, see Table 1. The table shows the target ID, the description of the target, the organism, the compounds and end points. In the case of the Beta-secretase 1, presents 2 157 compounds, and 3 360 endpoints. These separate data sets, (see Fig. 4), and the entire ChEMBLdb, are available either via ftp downloads or .xls files download, or via bespoke query interfaces, tailored to the requirements of the scientific communities with a specific interest in these research areas. 2.1.2. PDB Structures of -Secretases for DOCKING Studies The Protein Data Bank (PDB) www.pdb.org archive is the single worldwide repository of information about the 3D structures of large biological molecules, including proteins and nucleic acids, (see Fig. 5). These are the molecules of life that are found in all organisms including bacteria, yeast, plants, flies, other animals, and humans. Understanding the shape of a molecule helps to understand how it works. This knowledge can be used to help deduce a structure’s role in human health and disease, and in drug development. The structures in the archive range from tiny proteins and bits of DNA to complex molecular machines like the ribosome. 3D Fig. (3). ChEMBL database server. 12 Current Bioinformatics, 2011, Vol. 6, No. 1 Prado-Prado et al. drugs. The codes and activity for all compounds as well as the references used to collect them are depicted in SM1 of the supplementary material file. 4.2.2. ANN Models The ANN models are non-linear models useful to predict the biological activity of a large datasets of molecules. This technique is an alternative to linear methods such as LDA. (Fig. 9) depicts the networks maps for some of the ANN models. In general, at least one ANN of every types tested was statically significant. However, one must note that the Fig. (9). Topology of some ANN models trained in this work. Table 7. Comparison of Different ANNs Classification Models Model Train Stat. Validation profile active Non-active % Par. active Non-active % RBF 1255 849 59.65 Sn 628 425 59.64 1:1-326-1:1 2385 5915 71.27 Sp 1186 2964 71.42 68.92 Ac 69.04 PNN 34 2070 1.62 Sn 17 1036 1.61 15:15-10404-2-2:1 0 8300 100 Sp 0 4150 100 80.10 Ac 80.09 MLP 1812 292 86.12 Sn 909 144 86.32 15:15-11-1:1 1118 7182 86.53 Sp 556 3594 86.60 86.45 Ac 86.55 Linear 1926 178 91.54 Sn 973 80 92.40 15:15-1:1 745 7555 91.02 Sp 308 3842 92.58 91.13 Ac 92.54 Review of Bioinformatics and QSAR Studies Current Bioinformatics, 2011, Vol. 6, No. 1 13 profiles of each network indicate that these are highly nonlinear and complicated models. There are several different kinds of ANN and these include multilayer perceptron (MLP), radial basis functions (RBF) and PNNs; the latter ANN is a variant of RBF systems. In particular, PNN is a type of neural network that uses a kernel-based approximation to form an estimate of the probability density functions of classes in a classification problem [84]. 4.3. Results and Discussion The network found was LNN and it showed training performance higher than 91%. We compare different types of networks to obtain a better model; Table 7 shows the classification matrix of the different networks. Linear 15:15-1:1 was taken as the main network because it presented a wider range of variables, 15 inputs in the first layer and 15 neurons in second layer, and two sets of cases (Training and Validation). Another tested networks found were MLP 15:15-111:1, RBF 1:1-326-1:1 presented low accuracy and PNN 15:15-10404-2-2:1 had a very low percentage of non-active leading to possible errors in the model although its accuracy very good, see Table 7. In Fig. (10), we depict the ROCcurve [85] for LNN tested. Notably, almost model presented and an area under curve higher than 0.5 (the value for a random classifier). The vitality of this type of procedures developing ANN-QSAR models has been demonstrated befote [86]; see, for instance, the work of Fernandez and Caballero [87]. The same is true about the ANNs tested, where is illustrated ROC-curve of ANN LNN with an area higher than 0.93. To show how important is this result, we compared the present model with other model used to address the same problem. We processed our data with Artificial Neural Networks (ANNs) looking for a better model. In general, the ANN LNN tested was statically significant. 5. CONCLUSIONS Theoretical studies such as QSAR models have become a very useful tool in this context to substantially reduce time and resources consuming experiments. The functions of - secretase and its implication in Alzheimer's disease have triggered an active search for potent and selective -secretase inhibitors. In this paper we can see that the development of theoretical and QSAR models to study -secretase inhibitors are usually not many achieved so far, and most of these works present docking studies. Watching this situation we need to develop QSAR models with -secretase inhibitors. In this sense, QSAR could play an important role in studying these -secretase inhibitors. QSARs can be used as predictive tools for the development of molecules. In this work we developed a new ANN LNN model using the ModesLab descriptors, based on a large database using about 10,000 different drugs obtained from the ChEMBL server. ACKNOWLEDGEMENTS F. Prado-Prado thanks sponsorships for research position at the University of Santiago de Compostela from Angeles Alvariño, Xunta de Galicia. All authors acknowledge the Project 07CSA008203PR. REFERENCES [1] Querfurth HW, LaFerla FM. Alzheimer's disease. N Engl J Med 2010; 362: 329-44. [2] Goate A, Chartier-Harlin MC, Mullan M, et al. Segregation of a missense mutation in the amyloid precursor protein gene with familial Alzheimer's disease. Nature 1991; 349: 704-6. [3] Levy-Lahad E, Wasco W, Poorkaj P, et al. Candidate gene for the chromosome 1 familial Alzheimer's disease locus. Science 1995; 269: 973-7. [4] Sherrington R, Rogaev EI, Liang Y, et al. Cloning of a gene bearing missense mutations in early-onset familial Alzheimer's disease. Nature 1995; 375: 754-60. [5] Citron M, Westaway D, Xia W, et al. Mutant presenilins of Alzheimer's disease increase production of 42-residue amyloid betaprotein in both transfected cells and transgenic mice. Nat Med 1997; 3: 67-72. [6] Holcomb L, Gordon MN, McGowan E, et al. Accelerated Alzheimer-type phenotype in transgenic mice carrying both mutant Fig. (10). ROC Curve for classifier. 14 Current Bioinformatics, 2011, Vol. 6, No. 1 Prado-Prado et al. amyloid precursor protein and presenilin 1 transgenes. Nat Med 1998; 4: 97-100. [7] Cole SL, Vassar R. The Alzheimer's disease beta-secretase enzyme, BACE1. Mol Neurodegener 2007; 2: 22. [8] Cole SL, Vassar R. The basic biology of BACE1: a key therapeutic target for alzheimer's disease. Curr Genomics 2007; 8: 509-30. [9] Haass C, Schlossmacher MG, Hung AY, et al. Amyloid betapeptide is produced by cultured cells during normal metabolism. Nature 1992; 359: 322-5. [10] Seubert P, Oltersdorf T, Lee MG, et al. Secretion of beta-amyloid precursor protein cleaved at the amino terminus of the betaamyloid peptide. Nature 1993; 361: 260-3. [11] Fan Y, Unwalla R, Denny RA, et al. Insights for predicting bloodbrain barrier penetration of CNS targeted molecules using QSPR approaches. J Chem Inf Model 2010; 50: 1123-33. [12] Chou KC. Structural bioinformatics and its impact to biomedical science. Curr Med Chem 2004; 11: 2105-34. [13] Chou KC, Wei DQ, Du QS, Sirois S, Zhong WZ. Progress in computational approach to drug development against SARS. Curr Med Chem 2006; 13: 3263-70. [14] Todeschini R, Consonni V. Handbook of Molecular Descriptors. Wiley-VCH 2002. [15] Estrada E, Uriarte E. Recent advances on the role of topological indices in drug discovery research. Curr Med Chem 2001; 8: 157388. [16] González-Díaz H, González-Díaz Y, Santana L, Ubeira FM, Uriarte E. Proteomics, networks and connectivity indices. Proteomics 2008; 8: 750-78. [17] Nunez MB, Maguna FP, Okulik NB, Castro EA. QSAR modeling of the MAO inhibitory activity of xanthones derivatives. Bioorg Med Chem Lett 2004; 14: 5611-7. [18] Terada M, Inaba M, Yano Y, et al. Growth-inhibitory effect of a high glucose concentration on osteoblast-like cells. Bone 1998; 22: 17-23. [19] Santana L, Uriarte E, González-Díaz H, Zagotto G, Soto-Otero R, Mendez-Alvarez E. A QSAR model for in silico screening of MAO-A inhibitors. Prediction, synthesis, and biological assay of novel coumarins. J Med Chem 2006; 49: 1149-56. [20] Marrero-Ponce Y. Linear indices of the "molecular pseudograph's atom adjacency matrix": definition, significance-interpretation, and application to QSAR analysis of flavone derivatives as HIV-1 integrase inhibitors. J Chem Inf Comput Sci 2004; 44: 2010-26. [21] Vilar S, Santana L, Uriarte E. Probabilistic neural network model for the in silico evaluation of anti-HIV activity and mechanism of action. J Med Chem 2006; 49: 1118-24. [22] Marrero-Ponce Y, Khan MT, Casanola Martin GM, et al. Prediction of Tyrosinase Inhibition Activity Using Atom-Based Bilinear Indices. ChemMedChem 2007; 2: 449-78. [23] Casanola-Martin GM, Marrero-Ponce Y, Khan MT, et al. TOMOCOMD-CARDD descriptors-based virtual screening of tyrosinase inhibitors: evaluation of different classification model combinations using bond-based linear indices. Bioorg Med Chem 2007; 15: 1483-503. [24] Casanola-Martin GM, Marrero-Ponce Y, Khan MT, et al. Dragon method for finding novel tyrosinase inhibitors: Biosilico identification and experimental in vitro assays. Eur J Med Chem 2007; 42: 1370-81. [25] Caballero J, Fernandez M. Artificial neural networks from MATLAB in medicinal chemistry. Bayesian-regularized genetic neural networks (BRGNN): application to the prediction of the antagonistic activity against human platelet thrombin receptor (PAR-1). Curr Top Med Chem 2008; 8: 1580-605. [26] Duardo-Sanchez A, Patlewicz G, Lopez-Diaz A. Current topics on software use in medicinal chemistry: intellectual property, taxes, and regulatory issues. Curr Top Med Chem 2008; 8: 1666-75. [27] Gonzalez MP, Teran C, Saiz-Urra L, Teijeira M. Variable selection methods in QSAR: an overview. Curr Top Med Chem 2008; 8: 1606-27. [28] Gonzalez-Diaz H. Quantitative studies on Structure-Activity and Structure-Property Relationships (QSAR/QSPR). Curr Top Med Chem 2008; 8: 1554. [29] Gonzalez-Diaz H, Prado-Prado F, Ubeira FM. Predicting antimicrobial drugs and targets with the MARCH-INSIDE approach. Curr Top Med Chem 2008; 8: 1676-90. [30] Helguera AM, Combes RD, Gonzalez MP, Cordeiro MN. Applications of 2D descriptors in drug design: a DRAGON tale. Curr Top Med Chem 2008; 8: 1628-55. [31] Ivanciuc O. Weka machine learning for predicting the phospholipidosis inducing potential. Curr Top Med Chem 2008; 8: 1691-709. [32] Vilar S, Cozza G, Moro S. Medicinal chemistry and the molecular operating environment (MOE): application of QSAR and molecular docking to drug discovery. Curr Top Med Chem 2008; 8: 1555-72. [33] Wang JF, Wei DQ, Chou KC. Drug candidates from traditional chinese medicines. Curr Top Med Chem 2008; 8: 1656-65. [34] Wang JF, Wei DQ, Chou KC. Pharmacogenomics and personalized use of drugs. Curr Top Med Chem 2008; 8: 1573-9. [35] Torrens F, Castellano G. Topological charge-transfer indices: from small molecules to proteins. Curr Proteomics 2009: 204-13. [36] Concu R, Dea-Ayuela MA, Perez-Montoto LG, et al. 3D entropy and moments prediction of enzyme classes and experimentaltheoretic study of peptide fingerprints in Leishmania parasites. Biochimica et Biophysica Acta 2009; 1794: 1784-94. [37] Ivanciuc O. Machine learning Quantitative Structure-Activity Relationships (QSAR) for peptides binding to human amphiphysin-1 SH3 domain. Curr Proteomics 2009; 4: 289-302. [38] Vilar S, Gonzalez-Diaz H, Santana L, Uriarte E. A network-QSAR model for prediction of genetic-component biomarkers in human colorectal cancer. J Theor Biol 2009; 261: 449-58. [39] Giuliani A, Di Paola L, Setola R. Proteins as networks: a mesoscopic approach using haemoglobin molecule as case study. Curr Proteomics 2009; 6: 235-45. [40] Chou KC. Pseudo amino acid composition and its applications in bioinformatics, proteomics and system biology. Curr Proteomics 2009; 6: 262-74. [41] Chen J, Shen B. Computational analysis of amino acid mutation: a proteome wide perspective. Curr Proteomics 2009; 6: 228-34. [42] Zhong WZ, Zhan J, Kang P, Yamazaki S. Gender specific drug metabolism of PF-02341066 in rats--role of sulfoconjugation. Curr Drug Metab 2010; 11: 296-306. [43] Wang JF, Chou KC. Molecular modeling of cytochrome P450 and drug metabolism. Curr Drug Metab 2010; 11: 342-6. [44] Mrabet Y, Semmar N. Mathematical methods to analysis of topology, functional variability and evolution of metabolic systems based on different decomposition concepts. Curr Drug Metab 2010; 11: 315-41. [45] Martinez-Romero M, Vazquez-Naya JM, Rabunal JR, et al. Artificial intelligence techniques for colorectal cancer drug metabolism: ontology and complex network. Curr Drug Metab 2010; 11: 34768. [46] Khan MT. Predictions of the ADMET properties of candidate drug molecules utilizing different QSAR/QSPR modelling approaches. Curr Drug Metab 2010; 11: 285-95. [47] Gonzalez-Diaz H, Duardo-Sanchez A, Ubeira FM, et al. Review of MARCH-INSIDE & complex networks prediction of drugs: ADMET, anti-parasite activity, metabolizing enzymes and cardiotoxicity proteome biomarkers. Curr Drug Metab 2010; 11: 379-406. [48] Gonzalez-Diaz H. Network topological indices, drug metabolism, and distribution. Curr Drug Metab 2010; 11: 283-4. [49] Garcia I, Diop YF, Gomez G. QSAR & complex network study of the HMGR inhibitors structural diversity. Curr Drug Metab 2010; 11: 307-14. [50] Chou KC. Graphic rule for drug metabolism systems. Curr Drug Metab 2010; 11: 369-78. [51] Concu R, Podda G, Ubeira FM, Gonzalez-Diaz H. Review of QSAR models for enzyme classes of drug targets: Theoretical background and applications in parasites, hosts, and other organisms. Curr Pharm Design 2010; 16: 2710-23. [52] Estrada E, Molina E, Nodarse D, Uriarte E. Structural contributions of substrates to their binding to P-Glycoprotein. A TOPS-MODE approach. Curr Pharm Design 2010; 16: 2676-709. [53] Garcia I, Fall Y, Gomez G. QSAR, docking, and CoMFA studies of GSK3 inhibitors. Curr Pharm Design 2010; 16: 2666-75. [54] Gonzalez-Diaz H. QSAR and complex networks in pharmaceutical design, microbiology, parasitology, toxicology, cancer, and neurosciences. Curr Pharm Design 2010; 16: 2598-600. [55] Gonzalez-Diaz H, Romaris F, Duardo-Sanchez A, et al. Predicting drugs and proteins in parasite infections with topological indices of complex networks: theoretical backgrounds, applications, and legal issues. Curr Pharm Design 2010; 16: 2737-64. Review of Bioinformatics and QSAR Studies Current Bioinformatics, 2011, Vol. 6, No. 1 15 [56] Marrero-Ponce Y, Casanola-Martin GM, Khan MT, Torrens F, Rescigno A, Abad C. Ligand-based computer-aided discovery of tyrosinase inhibitors. Applications of the TOMOCOMD-CARDD method to the elucidation of new compounds. Curr Pharm Design 2010; 16: 2601-24. [57] Munteanu CR, Fernandez-Blanco E, Seoane JA, et al. Drug discovery and design for complex diseases through QSAR computational methods. Curr Pharm Design 2010; 16: 2640-55. [58] Roy K, Ghosh G. Exploring QSARs with Extended Topochemical Atom (ETA) indices for modeling chemical and drug toxicity. Curr Pharm Design 2010; 16: 2625-39. [59] Speck-Planche A, Scotti MT, de Paulo-Emerenciano V. Current pharmaceutical design of antituberculosis drugs: future perspectives. Curr Pharm Design 2010; 16: 2656-65. [60] Vazquez-Naya JM, Martinez-Romero M, Porto-Pazos AB, et al. Ontologies of drug discovery and design for neurology, cardiology and oncology. Curr Pharm Design 2010; 16: 2724-36. [61] Estrada E. On the topological sub-structural molecular design (TOSS-MODE) in QSPR/QSAR and drug design research. SAR QSAR Environ Res 2000; 11: 55-73. [62] Overington J. ChEMBL. An interview with John Overington, team leader, chemogenomics at the European Bioinformatics Institute Outstation of the European Molecular Biology Laboratory (EMBLEBI). Interview by Wendy A. Warr. J Comput Aided Mol Des 2009; 23: 195-8. [63] Burgos PV, Mardones GA, Rojas AL, et al. Sorting of the Alzheimer's disease amyloid precursor protein mediated by the AP-4 complex. Dev Cell 2010; 18: 425-36. [64] Gordon WR, Vardar-Ulu D, L'Heureux S, et al. Effects of S1 cleavage on the structure, surface export, and signaling activity of human Notch1 and Notch2. PLoS One 2009; 4: e6613. [65] Iserloh U, Wu Y, Cumming JN, et al. Potent pyrrolidineand piperidine-based BACE-1 inhibitors. Bioorg Med Chem Lett 2008; 18: 414-7. [66] Geschwindner S, Olsson LL, Albert JS, et al. Discovery of a novel warhead against beta-secretase through fragment-based lead generation. J Med Chem 2007; 50: 5903-11. [67] Kong GK, Adams JJ, Harris HH, et al. Structural studies of the Alzheimer's amyloid precursor protein copper-binding domain reveal how it binds copper ions. J Mol Biol 2007; 367: 148-61. [68] Shiba T, Kametaka S, Kawasaki M, et al. Insights into the phosphoregulation of beta-secretase sorting signal by the VHS domain of GGA1. Traffic 2004; 5: 437-48. [69] Poulsen SA, Watson AA, Fairlie DP, Craik DJ. Solution structures in aqueous SDS micelles of two amyloid beta peptides of A beta(128) mutated at the alpha-secretase cleavage site (K16E, K16F). J Struct Biol 2000; 130: 142-52. [70] Al-Nadaf A, Abu Sheikha G, Taha MO. Elaborate ligand-based pharmacophore exploration and QSAR analysis guide the synthesis of novel pyridinium-based potent beta-secretase inhibitory leads. Bioorg Med Chem 2010; 18: 3088-115. [71] Pandey A, Mungalpara J, Mohan CG. Comparative molecular field analysis and comparative molecular similarity indices analysis of hydroxyethylamine derivatives as selective human BACE-1 inhibitor. Mol Divers 2010; 14: 39-49. [72] Polgar T, Keseru GM. Virtual screening for beta-secretase (BACE1) inhibitors reveals the importance of protonation states at Asp32 and Asp228. J Med Chem 2005; 48: 3749-55. [73] Hetenyi C, Paragi G, Maran U, Timar Z, Karelson M, Penke B. Combination of a modified scoring function with two-dimensional descriptors for calculation of binding affinities of bulky, flexible ligands to proteins. J Am Chem Soc 2006; 128: 1233-9. [74] Moitessier N, Therrien E, Hanessian S. A method for induced-fit docking, scoring, and ranking of flexible ligands. Application to peptidic and pseudopeptidic beta-secretase (BACE 1) inhibitors. J Med Chem 2006; 49: 5885-94. [75] Jung HA, Oh SH, Choi JS. Molecular docking studies of phlorotannins from Eisenia bicyclis with BACE1 inhibitory activity. Bioorg Med Chem Lett 2010; 20: 3211-5. [76] Zuo Z, Luo X, Zhu W, et al. Molecular docking and 3D-QSAR studies on the binding mechanism of statine-based peptidomimetics with beta-secretase. Bioorg Med Chem 2005; 13: 2121-31. [77] Salmon SA, Watts JL. Minimum inhibitory concentration determinations for various antimicrobial agents against 1570 bacterial isolates from turkey poults. Avian Dis 2000; 44: 85-98. [78] Chou KC. Review: Structural bioinformatics and its impact to biomedical science. Curr Med Chem 2004; 11: 2105-34. [79] Chou KC, Wei DQ, Du QS, Sirois S, Zhong WZ. Review: progress in computational approach to drug development against SARS. Curr Med Chem 2006; 13: 3263-70. [80] Prado-Prado FJ, Gonzalez-Diaz H, de la Vega OM, Ubeira FM, Chou KC. Unified QSAR approach to antimicrobials. Part 3: first multi-tasking QSAR model for input-coded prediction, structural back-projection, and complex networks clustering of antiprotozoal compounds. Bioorg Med Chem 2008; 16: 5871-80. [81] Prado-Prado FJ, de la Vega OM, Uriarte E, Ubeira FM, Chou KC, Gonzalez-Diaz H. Unified QSAR approach to antimicrobials. 4. Multi-target QSAR modeling and comparative multi-distance study of the giant components of antiviral drug-drug complex networks. Bioorg Med Chem 2009; 17: 569-75. [82] Kubinyi H. Quantitative structure-activity relationships (QSAR) and molecular modelling in cancer research. J Cancer Res Clin Oncol 1990; 116: 529-37. [83] Prado-Prado FJ, Borges F, Perez-Montoto LG, Gonzalez-Diaz H. Multi-target spectral moment: QSAR for antifungal drugs vs. different fungi species. Eur J Med Chem 2009; 44(10): 4051-6. [84] Mosier PD, Jurs PC. QSAR/QSPR studies using probabilistic neural networks and generalized regression neural networks. J Chem Inf Comput Sci 2002; 42: 1460-70. [85] Lombardi G, Gramegna G, Cavanna C, Michelone G. Fluconazole vs amphotericin B: "in vitro" comparative evaluation of the minimal inhibitory concentration (MIC) against yeasts isolated from AIDS patients. Microbiologica 1990; 13: 201-6. [86] Prado-Prado FJ, Garcia-Mera X, Gonzalez-Diaz H. Multi-target spectral moment QSAR versus ANN for antiparasitic drugs against different parasite species. Bioorg Med Chem 2010; 18: 2225-31. [87] Fernandez M, Caballero J, Tundidor-Camba A. Linear and nonlinear QSAR study of N-hydroxy-2-[(phenylsulfonyl)amino] acetamide derivatives as matrix metalloproteinase inhibitors. Bioorg Med Chem 2006; 14: 4137-50. Received: October 10, 2010 Revised: October 20, 2010 Accepted: December 09, 2010 Current Bioinformatics 2011, 0, 000-000 1574-8936/11 $55.00+.00 © 2011 Bentham Science Review of Bioinformatics and Theoretical studies of Acetylcholinesterase inhibitors Manuel Escobar, Franco Fernández, Xerardo García-Mera and Francisco Prado-Prado* Department of Organic Chemistry, University of Santiago de Compostela, 15782, Spain. Abstract: Alzheimer’s disease is a complex disease, and no single “magic bullet” is likely to prevent or cure it. That’s why current treatments focus on several different aspects, including helping people maintain mental function; managing behavioral symptoms; and slowing, delaying, or preventing the disease. Four medications are approved by the U.S. Food and Drug Administration to treat Alzheimer’s. Donepezil, rivastigmine, and galantamine are used to treat mild to moderate Alzheimer’s. Memantine is used to treat moderate to severe Alzheimer’s. These drugs work by regulating neurotransmitters (the chemicals that transmit messages between neurons). Treatment of AD by ACh precursors and cholinergic agonists was ineffective or caused severe side effects. ACh hydrolysis by AChE causes termination of cholinergic neurotransmission. Therefore, compounds which inhibit AChE might significantly increase the levels of ACh depleted in AD. However, these drugs don’t change the underlying disease process and may help only for a few months to a few years. In this sense, quantitative structure-activity relationships (QSAR) could play an important role in studying these AChE inhibitors. QSAR models are necessary in order to guide the AChE synthesis. In this work, we revised different bioinformatics and theoretical studies of Acetylcholinesterase inhibitors, design and computational studies for a very large and heterogeneous series of AChE inhibitors. First, we review 2D QSAR, 3D QSAR, CoMFA, CoMSIA and Docking and new theoretical methodology with different compound to find out the structural requirements. Next, we revised QSAR studies using method of Linear Discriminant Analysis (LDA) in order to understand the essential structural requirement for binding with receptor for AChE inhibitors. Keywords: QSAR; CoMSIA; COMFA; topological indices; Molecular Docking; Acetylcholinesterase inhibitors; Alzheimer's disease (AD). INTRODUCTION Alzheimer's disease (AD) is a disorder that attacks the central nervous system through progressive degeneration of its neurons. AD occurs in around 10% of the elderly and, as yet, there is no known cure. Patients with this disease develop dementia which becomes more severe as the disease progresses [1]. It was suggested that symptoms of AD are caused by decrased of activity of cholinergicneocortical and hippocampal neurons. Treatment of AD by ACh precursors and cholinergic agonists was ineffective or caused severe side effects. ACh hydrolysis by AChE causes termination of cholinergic neurotransmission. Therefore, compounds which inhibit AChE might significantly increase the levels of ACh depleted in AD. Indeed, it was shown that AChE inhibitors improve the cognitive abilities of AD patients at early stages of the disease development [2-4]. Acetylcholinesterase (AChE) is key enzyme in the nervous system of animals. By rapid hydrolysis of the neurotransmitter, acetylcholine (ACh), AChE terminates neurotransmission at cholinergic synapses. It is a very fast enzyme, especially for a serine hydrolase, functioning at a rate approaching that of a diffusioncontrolled reaction. *Corresponding author: Prado-Prado, Francisco, Department of Organic Chemistry, Faculty of Pharmacy, University of Santiago de Compostela, 15782, Email: francisco.[email protected]. AChE inhibitors are among the key drugs approved by the FDA for management of Alzheimer's disease (AD) [5-8]. The powerful toxicity of organophosphorus (OP) poisons is attributed primarily to their potent AChE inhibitors. Solution of the three-dimensional (3D) structure of Torpedo californica acetylcholinesterase (TcAChE) in 1991 opened up new horizons in research on an enzyme that had already been the subject of intensive investigation [9]. The unanticipated structure of this extremely rapid enzyme, in which theactive site was found to be buried at the bottom of a deep and narrow gorge, lined by 14 aromatic residues (colored dark magenta), led to a revision of the views then held concerning substrate traffic, recognition and hydrolysis [10]. To understand how those aromatic residues behave with the enzyme, see Flexibility of aromatic residues in acetylcholinesterase . Solution of the 3D structure of acetylcholinesterase led to a series of theoretical and experimental studies, which took advantage of recent advances in theoretical techniques for treatment of proteins or enzimes, such as molecular dynamics and electrostatics and to sitedirected mutagenesis, utilizing suitable expression systems. Acetylcholinesterase hydrolysizes theneurotra nsmitter acetylcholine (ACh), producing choline and an acetate group. Acetylcholine directly binds Ser200(via its nucleophilic Oγ atom) within the catalytic triad (Ser200, His440, and Glu327) (ACh/TcAChE structure 2ace). The residues Trp84 and Phe330 are also Bioinformatics and Theoretical studies of AChE inhibitors Current Bioinformatics 2011, 0, 000-000 important in the ligand recognition. See also: AChE inhibitors and substrates , see F ig ur e 1 . Figure 1. Cholinergic Synapse: Key Enzyme in the Nervous System. In this part, Chemoinformatics and Bioinformatics methods may play an important role in the study of acetylcholinesterase inhibitors, Quantitative StructureActivity Relationships (QSAR) studies are used as predictive tools for the molecular development [11, 12]. Up to today, there are near 1600 molecular descriptors that, in principle, can be generalized and used to solve the former problem [13]. Many of these indices are known as molecular Topological Indices (TIs) or simply invariants of a molecular graph. Unfortunately, QSAR studies are generally based on databases considering only structurally parent compounds acting against one single microbial species. In a recent review, our group have discussed recent advances in the field [14]. In addition to QSAR, Bioinformatics and Chemoinformatics methods useful to study β-secretase may include techniques like Comparative Molecular Field Analysis (CoMFA), drugtarget Docking, Sequence Alignment (SA) or other methods. In a recent, preliminary review in the field published in Proteomics in 2008 was discussed the use of these methods but only from the point of view of proteins [15]. Almost all QSAR techniques are based on the use of molecular descriptors, which are numerical series that codify useful chemical information and enable correlations between statistical and biological properties [16, 17]. On the other hand, QSAR models can be used to explore the relationships between the structural spaces of compounds as inhibitors for specific enzymes, such as MAO inhibitors [18], HIV-1 integrase inhibitors [19], and/or protease inhibitors [20] or tyrosinase inhibitors [21-23]. In fact, recently, the field has moved from small molecules to proteins and other systems. For instance, González-Díaz et al. discussed the use of these methods but only from the point of view of proteins [15]. Later, some groups published different papers in one special issue on QSAR but also restricted to the field of protein and proteomics [24-30]. In other recent issue, guestedited by González-Díaz [31] appeared a series of papers devoted to QSAR/QSPR techniques for low-molecularweight drugs [31-40]. Most recently, Prado-Prado et al. [41] published a mt-QSAR for anti-parasitic drugs. This year was published other issue [42] focused on QSAR/QSPR models and graph theory used to approach Drug ADMET processes and Metabolomics [43-50]. Last, two of the most recent issues published have been focused on the discussions of the applications of QSAR in Pharmaceutical Design [51-60] and Bioinformatics [6170]. In the present work, we review previous works based on 2D-QSAR, 3D-QSAR, CoMFA, CoMSIA and Docking techniques, which studied different compounds to find out the structural requirements. Last, we carried out new QSAR studies using Linear Discriminant Analysis (LDA) method and the software ModesLab[71] in order to understand the essential structural requirement for binding with receptor for acetylcholinesterase inhibitors. The topics reviewed, discussed, and/or reported in this paper are: 1. Theoretical studies for acetylcholinesterase inhibitors. 1.1. Acaricidal and QSAR of monoterpenes against Tetranychus urticae. 1.2. Preparation, AChE activity, docking study, and 3D-QSAR. 1.3. Prediction of AChE inhibitors and characterization by machine learning methods. 1.4. 3D-QSAR Studies of Physostigmine Analogues as AChE Inhibitors. 1.5. 3D-QSAR Studies on Carbamates as ACeH Inhibitors. 1.6. Classification of drug considering their IC50 using linear programming. 1.7. An investigation of carbamates for AChE inhibition using 3D-QSAR. 1.8. Molecular docking and 3D-QSAR studies indanone as AChE inhibitors. 1.9. QSAR of tacrine derivatives against AChE activity using variable selections. 1.10. 3D QSAR studies of AChE inhibitors based on molecular docking. 1.11. Anchor-GRIND: Filling the Gap between Standard 3D QSAR-GRIND. 1.12. A Docking Score Function for Estimating Ligand-Protein Interactions. 1.13. Modulation of Binding Strength in Several Classes of Active Site Inhibitors. 1.14. Structure-based 3D QSAR and design of novel ACeH inhibitors. 2. New method for the study of new acetylcholinesterase inhibitors. 2.1. Preface to new QSAR study acetylcholinesterase inhibitors 2.2. Methods 2.3. Results and Discussion. 1. Theoretical studies for acetylcholinesterase inhibitors. In this section we updated the contents presented in our recent review published in Current Drugs Current Bioinformatics 2011, 0, 000-000 Prado-Prado, et al. Metabolim [72]. The high number of possible candidates to acetylcholinesterase inhibitors creates the necessity of Quantitative Structure-Activity Relationship models in order to guide the acetylcholinesterase inhibitor synthesis. In this work, we revised different computational studies for a very large and heterogeneous series of acetylcholinesterase. First, we review 3D QSAR, CoMFA, CoMSIA and Docking with different compound to find out the structural requirements for acetylcholinesterase inhibitors. Last, we carried out new QSAR studies using Linear Discriminant Analysis (LDA) method and the software ModesLab[71] in order to understand the essential structural requirement for binding with receptor for acetylcholinesterase inhibitors. 1.1. Acaricidal and QSAR of monoterpenes against Tetranychus urticae. Badawy et al. [73] reviewed the acaricidal activity of 12 monoterpenes against the two-spotted spider mite, Tetranychus urticae Koch, was examined using fumigation and direct contact application methods. Cuminaldehyde and (-)-linalool showed the highest fumigant toxicity with LC50 = 0.31 and 0.56 mg/l, respectively. The other monoterpenes exhibited a strong fumigant toxicity, the LC50 values ranging from 1.28 to 8.09 mg/l, except camphene, which was the least effective (LC50 = 61.45 mg/l). Based on contact activity, the results were rather different: menthol displayed the highest acaricidal activity (LC50 = 128.53 mg/l) followed by thymol (172.0 mg/l), geraniol (219.69 mg/l) and (-)-limonene (255.44 mg/l); 1-8-cineole, cuminaldehyde and (-)-linalool showed moderate toxicity. At 125 mg/l, (-)-Limonene and (-)-carvone caused the highest egg mortality among the tested compounds (70.6 and 66.9% mortality, respectively). In addition, the effect of molecular descriptors was also analyzed using the quantitative structure activity relationship (QSAR) procedure. The QSAR model showed excellent agreement between the estimated and experimentally measured toxicity parameter (LC50) for the tested monoterpenes and the fumigant activity increased significantly with the vapor pressure. The authors compared the results of the fumigant and contact toxicity assays of monoterpenes against T. urticae with the results of acetylcholinesterase (AChE) inhibitory effect revealed that some of the tested compounds showed a strong acaricidal activity and a potent AChE inhibitory activity, such as cuminaldehyde, (-)-linalool, (-)-limonene and menthol. However, other compounds such as (-)-carvone revealed a strong fumigant activity but a weak AChE inhibitory activity. 1.2. Preparation, AChE activity, docking study, and 3D-QSAR Synthesis and anticholinesterase activity of 4-aryl-4oxo-N-phenyl-2-aminylbutyramides, novel class of reversible, moderately potent cholinesterase inhibitors, are reported by Vitorovic-Todorovic et al.[74]. In this work the authors used a simple substituent variation on aroyl moiety changes anti-AChE activity for two orders of magnitude; also substitution and type of hetero ( ali) cycle in position 2 of butanoic moiety govern AChE/BChE selectivity. The most potent compounds showed mixed-type inhibition, indicating their binding to free enzyme and enzyme–substrate complex. Figure 2. Compounds 5(R), (a); 5(S), (b); 19(R), (c); and 19(S), (d), docked into the binding site of the AChE, highlightening the protein residues that form the main interactions with the different structural units of the inhibitors. Hydrogen bonds are represented as black dots. Residues are colored as follows: anionic site Trp 86 and Tyr 337, orange; PAS Tyr 72, Asp 74, Tyr 124, Ser 125 violet and Trp 286 magenta; acyl binding pocket Tyr 341, Phe 297 and Phe 338 green; catalytic triad residues His 447 and Ser 203 dark blue; ligand, turquoise. Alignment-independent 3D QSAR study on reported compounds, and compounds having similar potencies obtained from the literature, confirmed that alkyl substitution on aroyl moiety of molecules is requisite for inhibition activity. The presence of hydrophobic moiety at close distance from hydrogen bond acceptor has favorable influence on inhibition potency. Docking studies show that compounds probably bind in the middle of the AChE active site gorge, but are buried deeper inside BChE active site gorge, as a consequence of larger BChE gorge void, see Figure 2. Bioinformatics and Theoretical studies of AChE inhibitors Current Bioinformatics 2011, 0, 000-000 1.3. Prediction of AChE inhibitors and characterization by machine learning methods. Acetylcholinesterase (AChE) has become an important drug target and its inhibitors have proved useful in the symptomatic treatment of Alzheimer’s disease. This work Wei Lv et al. [75] explores several machine learning methods (support vector machine (SVM), k-nearest neighbor (k-NN), and C4.5 decision tree (C4.5 DT)) for predicting AChE inhibitors (AChEIs). A feature selection method is used for improving prediction accuracy and selecting molecular descriptors responsible for distinguishing AChEIs and non-AChEIs. The prediction accuracies are 76.3%w88.0% for AChEIs and 74.3%w79.6% for nonAChEIs based on the three kinds of machine learning methods. This work suggests that machine learning methods such as SVM are facilitating for predicting AChEIs potential of unknown sets of compounds and for exhibiting the molecular descriptors associated with AChEIs, see Figure 3. Figure 3. The structures of the misclassified AChEIs. 1.4. 3D-QSAR Studies of Physostigmine Analogues as AChE Inhibitors. Natural alkaloid Physostigmine is one of the most potent pseudo-irreversible inhibitor of Acetylcholinesterase. It was found to accelerate longterm memory process, but due to its short half life and variable bioavailability, has inconsistent clinical efficacy. Zaheer Ul-Haq et al. [76] proposed a 3DQSAR studies based on the comparative molecular field analysis and comparative molecular similarity indices analysis and were applied to a set of 40 Physostigmine derivatives which are divided into two classes: A and B. In this study was obtained a highly reliable and extensive dynamic QSAR model based on alignment procedure with co-crystallized Ganstigmine as template. The strategy yielded significant 3DQSAR models with the cross-validated q2 values 0.762 and 0.754 for comparative molecular field analysis and comparative molecular similarity indices analysis, respectively. Resulted models were validated by external set of eight compounds yielding high correlation coefficient r2 values of 0.730 and 0.720 for comparative molecular field analysis and comparative molecular similarity indices analysis, respectively. Furthermore, the analysis of comparative molecular field analysis and comparative molecular similarity indices analysis contour maps within the active site of AChE were conducted in order to understand the interactions between the receptor and the Physostigmine derivatives, see Figure 4. Current Bioinformatics 2011, 0, 000-000 Prado-Prado, et al. Figure 4. 3D CoMSIA steric (A), CoMSIA electrostatic (B), CoMSIA hydrophobic (C), CoMSIA H-bond donor and acceptor (D) contour maps. This class of study will facilitate the rational design of more potent Physostigmine compounds which might have better activity and reduce toxicity for the treatment of Alzheimer disease. 1.5. 3D-QSAR Studies on Carbamates as ACeH Inhibitors In view of the nonavailability of complete X-ray structure of carbamates cocrystallized with AChE enzyme, the 3D-QSAR model development based on cocrystallized conformer (CCBA) as well as docked conformerbased alignment (DCBA) is not feasible. Therefore, the only two alternatives viz. pharmacophore and maximum common substructure-based alignments are left for the 3D-QSAR comparative molecular field analyses (CoMFA) and comparative molecular similarity indices analyses (CoMSIA) model development. Shailendra S. Chaudhaery et al. [77] in the presents iin this study, a 3D-QSAR models that have been developed using both alignment methods, where CoMFA and CoMSIA models based on pharmacophore-based alignment were in good agreement with each other and demonstrated significant superiority over MCS-based alignment in terms of leave-one-out (LOO) cross-validated q2 values of 0.573 and 0.723 and the r2 values of 0.972 and 0.950, respectively. The validation of the best CoMFA and CoMSIA models based on pharmacophore (Hip-Hop)- based alignment on a test set of 17 compounds provided significant predictive r2 [r2 pred(test)] of 0.614 and 0.788, respectively. The authors presented a contour map analyses, that revealed the relative importance of steric, electrostatic, and hydrophobicity for AChE inhibition activity, see Figure 5. However, hydrophobic factor plays a major contribution to the AChE inhibitory activity modulation which is in strong agreement with the fact that the AChE is having a wide active site gorge (∼20 Å) occupied by a large number of hydrophobic amino acid residues. Figure 5. The docked bioactive conformation of the most active compound 1 into the active site of the AChE enzyme generated by GOLD (version 3.0.1). Representation of hydrophobic (A) and electrostatic (B) zones around the most active compound 1. 1.6. Classification of drug considering their IC50 using linear programming A priori analysis of the activity of drugs on the target protein by computational approaches can be useful in narrowing down drug candidates for further experimental tests. Currently, there are a large number of computational methods that predict the activity of drugs on proteins. In this study, Pelin Armutlu et al. [78] approach the activity prediction problem as a classification problem and, they aim to improve the classification accuracy by introducing an algorithm that combines partial least squares regression with mixedinteger programming based hyper-boxes classification method, where drug molecules are classified as low active or high active regarding their binding activity (IC50 values) on target proteins. In this work they also aim to determine the most significant molecular Bioinformatics and Theoretical studies of AChE inhibitors Current Bioinformatics 2011, 0, 000-000 descriptors for the drug molecules, see Figure 6. First analyzing the activities of widely known inhibitor datasets including Acetylcholinesterase (ACHE), Benzodiazepine Receptor (BZR), Dihydrofolate Reductase (DHFR), Cyclooxygenase-2 (COX-2) with known IC50 values. The results at this stage proved that their approach consistently gives better classification accuracies compared to 63 other reported classification methods such as SVM, Naïve Bayes, where they were able to predict the experimentally determined IC50 values with a worst case accuracy of 96%. To further test applicability of this approach first the authors created dataset for Cytochrome P450 C17 inhibitors and then predicted their activities with 100% accuracy. The authors concluded results indicate that this approach can be utilized to predict the inhibitory effects of inhibitors based on their molecular descriptors. Figure 6. Outline of classification approach. 1.7. An investigation of carbamates for AChE inhibition using 3D-QSAR analysis Kuldeep K. Roy et al. in order to identify the essential structural features and physicochemical properties for AChE inhibitory activity in some carbamate derivatives, the systematic QSAR studies (CoMFA, advance CoMFA and CoMSIA) have been carried out on a series of (total 78 molecules) taking 52 and 26 molecules in training and test set, respectively. Statistically significant 3D-QSAR (three-dimensional Quantitative Structure Activity Relationship) models were developed on training setmolecules using CoMFA and CoMSIA and validated against test set compounds. The highly predictive models (CoMFA q2 = 0.733, r2 = 0.967, predictive r2 = 0.732, CoMSIA q2 = 0.641, r2 = 0.936, predictive r2 = 0.812) well explained the variance in binding affinities both for the training and the test set compounds. The generated models suggest that steric, electrostatic and hydrophobic interactions play an important role in describing the variation in binding affinity. In particular the carbamoyl nitrogen should be more electropositive; substitutions on this nitrogen should have high steric bulk and hydrophobicity while the amino nitrogen should be electronegative in order to have better activity, see Figure 7. These studies may provide important insights into structural variations leading to the development of novel AChE inhibitorswhich may be useful in the development of novelmolecules for the treatment of Alzheimer’s disease. Figure 7. Contour map of the best CoMFA model using Tripos standard field and contour map of the best CoMSIA model using steric, electrostatic and hydrophobic fields. 1.8. Molecular docking and 3D-QSAR studies indanone as AChE inhibitors Liang-liang Shen et al.[79] explore the binding mode of 2-substituted 1-indanone derivatives with acetylcholinesterase (AChE) and provide hints for the future design of new derivatives with higher potency and specificity. The GOLD-docking conformations of the compounds in the active site of the enzyme were used in subsequent studies were the method used by the autors. The highly reliable and predictive 3D-QSAR models were achieved by comparative molecular field analysis (CoMFA) and comparative molecular similarity analysis (CoMSIA) methods. The predictive capabilities of the models were validated by an external test set. Moreover, the stabilities of the 3D-QSAR models were verified by the leave-4-out crossvalidation method. The CoMFA and CoMSIA models were constructed successfully with a good crossvalidated coefficient (q2) and a non-cross-validated coefficient (r2). The q2 and r2 obtained from the leave-1out cross validation method were 0.784 and 0.974 in the Current Bioinformatics 2011, 0, 000-000 Prado-Prado, et al. CoMFA model and 0.736 and 0.947 in the CoMSIA model, respectively. The coefficient isocontour maps obtained from these models were compatible with the geometrical and physicochemical properties of AChE see Figure 8. A conclusion of this method was that the contour map demonstrated the binding affinity could be enhanced when the small protonated nitrogen moiety was replaced by a more hydrophobic and bulky group with a highly partial positive charge. The authors in this study provide a better understanding of the interaction between the inhibitors and AChE, which is helpful for the discovery of new compounds with more potency and selective activity. Figure 8. Alignment of all 45 compounds based on GOLD results. This image was generated with the Base program in SYBYL version 7.0. 1.9. QSAR of tacrine derivatives against AChE activity using variable selections Jung et al.[80] developed a diverse approach to the quantitative structure–activity relationship (QSAR) of tacrine derivatives against acetylcholinesterase (AChE) activity and was studied using variable selections of stepwise multiple linear regression (MLR), genetic algorithm (GA)- MLR, and simulated annealing (SA)- MLR. AChE activity (logRA) of tacrine derivatives was expressed with acceptable explanation (95.5–95.9%) and good predictive power (94.5–95.2%), respectively, in the models, see Table 1. The best equation was obtained from simulated annealing (SA) MLR with greater explanatory capability and better prediction, with a smaller standard error than other methods. The resulting models with the given descriptors illustrate the significant roles of hydrophobic and electrostatic interaction on increasing AChE activity, but hydrophilic and topological feature of molecules were shown to decrease AChE activity. Table 1. AChE activity of tacrine derivatives. Mol. ID Obs. Pred. Pred. Pred. Pred. Resid. Resid. 1 -1.37 -1.10 -0.99 -1.24 -0.17 0.14 -0.09 2 -1.41 -1.12 -1.48 -0.78 -0.27 -0.37 -0.58 3 -1.41 -1.40 -1.45 -1.24 -0.29 0.06 -0.17 4 -1.23 -1.40 -1.45 -1.19 -0.01 0.04 -0.21 5 -1.19 -1.26 -1.49 -1.19 0.17 0.22 -0.03 6 -1.26 -1.27 -1.44 -1.10 0.07 0.29 -0.08 7 -1.33 -1.15 -1.48 -1.10 0.01 0.17 -0.15 8 -0.71 -1.33 -1.37 -1.24 -0.17 0.14 -0.09 9b -0.43 -1.36 -1.37 -0.99 0.62 0.65 0.28 10b -0.68 -1.18 -1.15 -1.16 0.93 0.94 0.72 11b -0.80 -1.39 -1.51 -0.95 0.50 0.47 0.27 12 -0.85 -1.18 -1.08 -1.19 0.58 0.70 0.39 13 -0.74 -0.65 -0.58 -1.31 0.33 0.23 0.46 14b -0.74 -1.31 -1.19 -0.88 -0.09 -0.16 0.13 1.10. 3D QSAR studies of AChE inhibitors based on molecular docking and CoMFA Akula et al. [81] were performed a Threedimensional quantitative structure–activity relationship (3D QSAR) studies on acetylcholinesterase (AChE) inhibitors, based on molecular docking scores obtained by using FlexX and FlexiDock and comparative molecular field analysis (CoMFA). The docking scores were used as molecular descriptors along with the steric and electrostatic field values of CoMFA, for partial least square (PLS) analysis. The high leave one out (LOO) cross-validated correlation coefficient (q2 = 0.714) reveals that the model is a useful tool for the prediction of test set as well as newly designed structures against AChE activity. The superimposed CoMFA models on the receptor site of AChE are guiding the design of potential inhibitory structures directed against AChE activity, see Figure 9. Bioinformatics and Theoretical studies of AChE inhibitors Current Bioinformatics 2011, 0, 000-000 [51] Nassif H, Al-Ali H, Khuri S, Keirouz W. Prediction of protein-glucose binding sites using support vector machines. Proteins 2009; 77(1): 121-32. [52] Del Rio A, Baldi BF, Rastelli G. Activity prediction and structural insights of extracellular signal-regulated kinase 2 inhibitors with molecular dynamics simulations. Chemical biology & drug design 2009; 74(6): 630-5. [53] Garcia I, Fall Y, Gomez G. QSAR, docking, and CoMFA studies of GSK3 inhibitors. Curr Pharm Des 2010; 16(24): 2666-75. [54] Le T, Tseng TT, Saier MH, Jr. Flexible programs for the prediction of average amphipathicity of multiply aligned homologous proteins: application to integral membrane transport proteins. Mol Membr Biol 1999; 16(2): 173-9. [55] Gonzalez-Diaz H, Romaris F, Duardo-Sanchez A, Perez-Mototo LG, Prado-Prado F, Patlewicz G, Ubeira FM. Predicting drugs and proteins in parasite infections with topological indices of complex networks: theoretical backgrounds, aplications, and legal issues. Curr Pharm Des 2010; 16(24): 2737-64. [56] Marrero-Ponce Y, Casanola-Martin GM, Khan MT, Torrens F, Rescigno A, Abad C. Ligand-Based Computer-Aided Discovery of Tyrosinase Inhibitors. Applications of the TOMOCOMD-CARDD Method to the Elucidation of New Compounds. Curr Pharm Des 2010; 16(24): 2601-24. [57] Munteanu CR, Fernandez-Blanco E, Seoane JA, Izquierdo-Novo P, Rodriguez-Fernandez JA, PrietoGonzalez JM, Rabunal JR, Pazos A. Drug discovery and design for complex diseases through QSAR computational methods. Curr Pharm Des 2010; 16(24): 2640-55. [58] Roy K, Ghosh G. Exploring QSARs with Extended Topochemical Atom (ETA) Indices for Modeling Chemical and Drug Toxicity. Curr Pharm Des 2010; 16(24): 2625-39. [59] Speck-Planche A, Scotti MT, de PauloEmerenciano V. Current pharmaceutical design of antituberculosis drugs: future perspectives. Curr Pharm Des 2010; 16(24): 2656-65. [60] Veluraja K, Seethalakshmi AN. Dynamics of sialyl Lewis(a) in aqueous solution and prediction of the structure of the sialyl Lewis(a)-SelectinE complex. J Theor Biol 2008; 252(1): 15-23. [61] Bhattacharjee B, Jayadeepa RM, Banerjee S, Joshi J, Middha SK, Mole JP, Samuel J. Review of Complex Network and Gene Ontology in pharmacology approaches: Mapping natural compounds on potential drug target Colon Cancer network. Current Bioinformatics 2011; 6(1): 44-52. [62] Chiş O, Dumitru O, Concu R, Shen B. Reviewing Yeast Network and report of new Stochastic-Credibility cell cycle models. Current Bioinformatics 2011; 6(1): 35-43. [63] Dave K, Banerjee A. Bioinformatics analysis of functional relations between CNPs regions. Current Bioinformatics 2011; 6(1): 122-8. [64] Duardo-Sanchez A, Patlewicz G, González-Díaz H. A Review of Network Topological Indices from Chem-Bioinformatics to Legal Sciences and back. Current Bioinformatics 2011; 6(11): 53-70. [65] García I, Fall Y, Gómez G. Trends in Bioinformatics and Chemoinformatics of Vitamin D analogues and their protein targets. Current Bioinformatics 2011; 6(1): 16-24. [66] Ivanciuc T, Ivanciuc O, Klein DJ. Network-QSAR with Reaction Poset Quantitative SuperstructureActivity Relationships (QSSAR) for PCB Chromatographic Properties. Current Bioinformatics 2011; 6(1): 25-34. [67] Prado-Prado F, Escobar-Cubiella M, GarcíaMera X. Review of Bioinformatics and QSAR studies of β-secretase inhibitors. Current Bioinformatics 2011; 6(1): 3-15. [68] Riera-Fernández P, Munteanu CR, PedreiraSouto N, Duardo-Sanchez A, González-Díaz H. Definition of Markov-Harary Invariants and Review of Classic Topological Indices and Databases of Biology, Parasitology, Technology, and Social-Legal Networks. Current Bioinformatics 2011; 6(1): 94-121. [69] Speck-Planche A, Cordeiro MNDS. Application of Bioinformatics for the search of novel anti-viral therapies: Rational design of anti-herpes agents. Current Bioinformatics 2011; 6(1): 81-93. [70] Wan SB, Hu LL, Niu S, Wang K, Cai YD, Lu WC, Chou KC. Identification of multiple subcellular locations for proteins in budding yeast. Current Bioinformatics 2011; 6(1): 71-80. [71] Estrada E. On the topological sub-structural molecular design (TOSS-MODE) in QSPR/QSAR and drug design research. SAR QSAR Environ Res 2000; 11(1): 55-73. [72] Brady CP, Dowd AJ, Tort J, Roche L, Condon B, O'Neill SM, Brindley PJ, Dalton JP. The cathepsin L-like proteinases of liver fluke and blood fluke parasites of the trematode genera Fasciola and Schistosoma. Biochem Soc Trans 1999; 27(4): 740-5. [73] Badawy ME, El-Arami SA, Abdelgaleil SA. Acaricidal and quantitative structure activity relationship of monoterpenes against the two-spotted spider mite, Tetranychus urticae. Exp Appl Acarol 2010; 52(3): 261-74. [74] Vitorovic-Todorovic MD, Juranic IO, Mandic LM, Drakulic BJ. 4-Aryl-4-oxo-N-phenyl-2aminylbutyramides as acetyland butyrylcholinesterase inhibitors. Preparation, anticholinesterase activity, docking study, and 3D structure-activity relationship based on molecular interaction fields. Bioorganic & medicinal chemistry 2010; 18(3): 1181-93. Current Bioinformatics 2011, 0, 000-000 Prado-Prado, et al. [75] Lv W, Xue Y. Prediction of acetylcholinesterase inhibitors and characterization of correlative molecular descriptors by machine learning methods. Eur J Med Chem 2010; 45(3): 1167-72. [76] Ul-Haq Z, Mahmood U, Jehangir B. Ligandbased 3D-QSAR studies of physostigmine analogues as acetylcholinesterase inhibitors. Chemical biology & drug design 2009; 74(6): 571-81. [77] Chaudhaery SS, Roy KK, Saxena AK. Consensus superiority of the pharmacophore-based alignment, over maximum common substructure (MCS): 3D-QSAR studies on carbamates as acetylcholinesterase inhibitors. Journal of chemical information and modeling 2009; 49(6): 1590-601. [78] Armutlu P, Ozdemir ME, Uney-Yuksektepe F, Kavakli IH, Turkay M. Classification of drug molecules considering their IC50 values using mixed-integer linear programming based hyper-boxes method. BMC Bioinformatics 2008; 9: 411. [79] Shen LL, Liu GX, Tang Y. Molecular docking and 3D-QSAR studies of 2-substituted 1-indanone derivatives as acetylcholinesterase inhibitors. Acta pharmacologica Sinica 2007; 28(12): 2053-63. [80] Ahn JH, Shin MS, Jung SH, Kim JA, Kim HM, Kim SH, Kang SK, Kim KR, Rhee SD, Park SD, Lee JM, Lee JH, Cheon HG, Kim SS. Synthesis and structure-activity relationship of novel indene N-oxide derivatives as potent peroxisome proliferator activated receptor gamma (PPARgamma) agonists. Bioorg Med Chem Lett 2007; 17(18): 5239-44. [81] Akula N, Lecanu L, Greeson J, Papadopoulos V. 3D QSAR studies of AChE inhibitors based on molecular docking scores and CoMFA. Bioorganic & medicinal chemistry letters 2006; 16(24): 6277-80. [82] Fontaine F, Pastor M, Zamora I, Sanz F. Anchor-GRIND: filling the gap between standard 3D QSAR and the GRid-INdependent descriptors. J Med Chem 2005; 48(7): 2687-94. [83] Guo J, Hurley MM, Wright JB, Lushington GH. A docking score function for estimating ligand-protein interactions: application to acetylcholinesterase inhibition. J Med Chem 2004; 47(22): 5492-500. [84] Martin-Santamaria S, Munoz-Muriedas J, Luque FJ, Gago F. Modulation of binding strength in several classes of active site inhibitors of acetylcholinesterase studied by comparative binding energy analysis. J Med Chem 2004; 47(18): 4471-82. [85] Sippl W, Contreras JM, Parrot I, Rival YM, Wermuth CG. Structure-based 3D QSAR and design of novel acetylcholinesterase inhibitors. J Comput Aided Mol Des 2001; 15(5): 395-410. [86] Salmon SA, Watts JL. Minimum inhibitory concentration determinations for various antimicrobial agents against 1570 bacterial isolates from turkey poults. Avian Dis 2000; 44(1): 85-98. [87] Chou KC. Review: Structural bioinformatics and its impact to biomedical science. Curr Med Chem 2004; 11: 2105-34. [88] Chou KC, D. Q. Wei, Q. S. Du, S. Sirois, and W. Z. Zhong. . Review: Progress in computational approach to drug development against SARS. Curr Med Chem 2006; 13: 3263-70. [89] Prado-Prado FJ, Gonzalez-Diaz H, de la Vega OM, Ubeira FM, Chou KC. Unified QSAR approach to antimicrobials. Part 3: first multi-tasking QSAR model for input-coded prediction, structural back-projection, and complex networks clustering of antiprotozoal compounds. Bioorg Med Chem 2008; 16(11): 5871-80. [90] Prado-Prado FJ, de la Vega OM, Uriarte E, Ubeira FM, Chou KC, Gonzalez-Diaz H. Unified QSAR approach to antimicrobials. 4. Multi-target QSAR modeling and comparative multi-distance study of the giant components of antiviral drug-drug complex networks. Bioorg Med Chem 2009; 17: 569–75. [91] Kubinyi H. Quantitative structure-activity relationships (QSAR) and molecular modelling in cancer research. J Cancer Res Clin Oncol 1990; 116(6): 529-37. [92] Prado-Prado FJ, Borges F, Perez-Montoto LG, Gonzalez-Diaz H. Multi-target spectral moment: QSAR for antifungal drugs vs. different fungi species. Eur J Med Chem 2009. [93] Overington J. ChEMBL. An interview with John Overington, team leader, chemogenomics at the European Bioinformatics Institute Outstation of the European Molecular Biology Laboratory (EMBL-EBI). Interview by Wendy A. Warr. J Comput Aided Mol Des 2009; 23(4): 195-8. [94] Estrada E, Patlewicz G, Chamberlain M, Basketter D, Larbey S. Computer-aided knowledge generation for understanding skin sensitization mechanisms: the TOPS-MODE approach. Chem Res Toxicol 2003; 16(10): 1226-35. [95] Estrada E, Uriarte E, Gutierrez Y, Gonzalez H. Quantitative structure-toxicity relationships using TOPS-MODE. 3. Structural factors influencing the permeability of commercial solvents through living human skin. SAR QSAR Environ Res 2003; 14(2): 145-63. [96] Estrada E, Gonzalez H. What are the limits of applicability for graph theoretic descriptors in QSPR/QSAR? Modeling dipole moments of aromatic compounds with TOPS-MODE descriptors. J Chem Inf Comput Sci 2003; 43(1): 75-84. [97] Estrada E, Vilar S, Uriarte E, Gutierrez Y. In silico studies toward the discovery of new anti-HIV nucleoside compounds with the use of TOPS-MODE and 2D/3D connectivity indices. 1. Pyrimidyl derivatives. J Chem Inf Comput Sci 2002; 42(5): 1194203. [98] Estrada E. Quantum-chemical foundations of the topological substructural molecular design. J Phys Chem A 2008; 112(23): 5208-17. Bioinformatics and Theoretical studies of AChE inhibitors Current Bioinformatics 2011, 0, 000-000 [99] Yoshii F, Hirono S. Construction of a quantitative three-dimensional model for odor quality using comparative molecular field analysis (CoMFA). Chem Senses 1996; 21(2): 201-10. [100] Gonzalez-Diaz H, Saiz-Urra L, Molina R, Santana L, Uriarte E. A model for the recognition of protein kinases based on the entropy of 3D van der Waals interactions. Journal of proteome research 2007; 6(2): 904-8. [101] Alvarez-Ginarte YM, Marrero-Ponce Y, RuizGarcia JA, Montero-Cabrera LA, Vega JM, Noheda Marin P, Crespo-Otero R, Zaragoza FT, GarciaDomenech R. Applying pattern recognition methods plus quantum and physico-chemical molecular descriptors to analyze the anabolic activity of structurally diverse steroids. J Comput Chem 2007. [102] Morales AH, Rodríguez-Borges JE, García-Mera X, Fernández F, Dias-Sueiro-Cordeiro MN. Probing the Anticancer Activity of Nucleoside Analogues: A QSAR Model Approach Using an Internally Consistent Training Set. J Med Chem 2007; 50: 1537-45. Original article 2D MI-DRAGON: A new predictor for proteineligands interactions and theoreticexperimental studies of US FDA drug-target network, oxoisoaporphine inhibitors for MAO-A and human parasite proteins Francisco Prado-Prado a , * , Xerardo García-Mera a , Manuel Escobar a , Eduardo Sobarzo-Sánchez b , Matilde Yañez c , Pablo Riera-Fernandez d , Humberto González-Díaz d a Department of Organic Chemistry, Faculty of Pharmacy, University of Santiago de Compostela (USC) 15782, Spain b Department of Pharmaceutical Technology, Faculty of Pharmacy, USC 5782, Spain c Department of Pharmacology, Faculty of Pharmacy, USC 15782, Spain d Department of Microbiology and Parasitology, Faculty of Pharmacy, USC 15782, Spain article info Article history: Received 20 June 2011 Received in revised form 22 September 2011 Accepted 26 September 2011 Available online 1 October 2011 Keywords: DrugeProtein interaction complex networks Protein structure networks Multi-target QSAR Markov model MAO A inhibitors abstract There are many pairs of possible DrugeProteins Interactions that may take place or not (DPIs/nDPIs) between drugs with high affinity/non-affinity for different proteins. This fact makes expensive in terms of time and resources, for instance, the determination of all possible ligandseprotein interactions for a single drug. In this sense, we can use Quantitative StructureeActivity Relationships (QSAR) models to carry out rational DPIs prediction. Unfortunately, almost all QSAR models predict activity against only one target. To solve this problem we can develop multi-target QSAR (mt-QSAR) models. In this work, we introduce the technique 2D MI-DRAGON a new predictor for DPIs based on two different well-known software. We use the software MARCH-INSIDE (MI) to calculate 3D structural parameters for targets and the software DRAGON was used to calculated 2D molecular descriptors all drugs showing known DPIs present in the Drug Bank (US FDA benchmark dataset). Both classes of parameters were used as input of different Artificial Neural Network (ANN) algorithms to seek an accurate non-linear mt-QSAR predictor. The best ANN model found is a Multi-Layer Perceptron (MLP) with profile MLP 21:21-31-1:1. This MLP classifies correctly 303 out of 339 DPIs (Sensitivity ¼89.38%) and 480 out of 510 nDPIs (Specificity ¼94.12%), corresponding to training Accuracy ¼92.23%. The validation of the model was carried out by means of external predicting series with Sensitivity ¼92.18% (625/678 DPIs; Specificity ¼90.12% (730/780 nDPIs) and Accuracy ¼91.06%. 2D MI-DRAGON offers a good opportunity for fast-track calculation of all possible DPIs of one drug enabling us to re-construct large drug-target or DPIs Complex Networks (CNs). For instance, we reconstructed the CN of the US FDA benchmark dataset with 855 nodes 519 drugs þ336 targets). We predicted CN with similar topology (observed and predicted values of average distance are equal to 6.7 vs. 6.6). These CNs can be used to explore large DPIs databases in order to discover both new drugs and/or targets. Finally, we illustrated in one theoreticexperimental study the practical use of 2D MI-DRAGON. We reported the prediction, synthesis, and pharmacological assay of 10 different oxoisoaporphines with MAO-A inhibitory activity. The more active compound OXO5 presented IC 50 ¼0.00083 m M, notably better than the control drug Clorgyline. 2011 Elsevier Masson SAS. All rights reserved. 1. Introduction The prediction of interactions between organic compounds to form drug-target pairs (DPIs) is a source piece on the combination of bioinformatics and proteome research towards drug discovery. Therefore, there is a strong incentive to develop new methods capable of detecting these potential drug-target interactions efficiently [1]. In this sense, we can use Quantitative StructureeActivity Relationships (QSAR) models [2] to carry out rational DPIs prediction. However, the use of QSAR techniques to predict DPIs is an area less explored until now. Classic QSAR models are equations that connect the structure of the drug, expressed by means of molecular descriptors, with the biological function. In classic QSAR studies we can use different free or commercially *Corresponding author. Faculty of Pharmacy, University of Santiago de Compostela, 15782, Spain. Fax: þ34 981 594912. E-mail address: francisco.pr[email protected] (F. Prado-Prado). Contents lists available at SciVerse ScienceDirect European Journal of Medicinal Chemistry journal homepage: http://www.elsevier.com/locate/ejmech 0223-5234/$ esee front matter 2011 Elsevier Masson SAS. All rights reserved. doi:10.1016/j.ejmech.2011.09.045 European Journal of Medicinal Chemistry 46 (2011) 5838e5851 available software to calculate structural parameters of the drug. Some of the more known software we can use to reach this goal are: DRAGON, CODESSA [3], MODES-LAB [4], TOMO-COMD [5], and MARCH-INSIDE (MI) [6]. The software DRAGON is one of the more complete calculating more than 1600 descriptors for drug structure including as zero- (0D) one- (1D), two- (2D), three-dimensional (3D) parameters. Unfortunately almost all QSAR models are able to predict the activity of drugs against only one target. To solve this problem we can develop multitarget QSAR (mt-QSAR) models to predict DPIs [7]. One way to develop this class mt-QSAR is incorporating into the QSAR equation parameters of the structure of the target (protein, DNA, RNA, etc.) in addition to the structural parameters of the drug present in classic QSAR. Ideally, we should use structural parameters of both drug and target calculated with the same software in order to standardize the entire algorithm. Disappointedly, many of known software used for QSAR have been developed to calculate only molecular indices of drugs. In any case, the theory used in these software packages can be easily extended to calculate molecular descriptors of targets. Consequently, it is only a matter of time that the authors of these programs develop new upgrades incorporating also molecular descriptors of the targets. That is the case of TOMO-COMD and MI; both have incorporated the calculation of drug and target structural parameters. In any case, until the best of our knowledge only MI has been used to develop and publish Fig. 1. Flowchart of all steps given in this work to develop the new model. F. Prado-Prado et al. / European Journal of Medicinal Chemistry 46 (2011) 5838e5851 5839 a QSAR predictor for DPIs. For instance, very recently we developed MIND-BEST [8] and NL MIND-BEST [9] two predictors based on structural parameters of drugs and proteins calculated with software MI. As was mentioned in the previous paragraph we can seek a QSAR predictor for DPIs using molecular descriptors of both drug and target. Nevertheless, many of the aforementioned software can calculate only structural parameters for drugs and only a few of them have been upgraded to incorporate parameters of targets. In this sense, we propose to seek QSAR models to predict DPIs using two different software packages (one for drug parameters and other for target ones) as an alternative to the yet limited use of one program only. In a last step, we can use this new type of QSAR model for the prediction of all possible DPIs/nDPIs (non-drugprotein interaction) in the global set of relationships between protein targets and all drugs to form the complex network (CN) of all DPIs. This type of CN for DPIs may become of the major relevance for the discovery of new drugs and targets as well. For instance, Yildirim, et al. [10] have built a bipartite graph composed of US Food and Drug Administration (US FDA) approved drugs and proteins linked by drug-target binary associations. The resulting network connects most drugs into a highly interlinked giant component, with strong local clustering of drugs of similar types according to Anatomical Therapeutic Chemical classification. Topological analyses of this network quantitatively showed an overabundance of ’follow-on’drugs, that is, drugs that target already targeted proteins. In this work, we developed for first time one mt-QSAR predictor of DPIs called 2D MI-DRAGON. The new technique combines the software DRAGON [11] to calculate structural parameters of drugs with software MI to calculate structural parameters of the target. The technique 2D MI-DRAGON uses these descriptors as input of one Artificial Neural Network (ANN) to seek the model. Both training and validation of the model was carried out by means of learning and external predicting series containing structural parameters for all DPIs present in the Drug Bank (US FDA benchmark dataset) downloaded from Drug Bank [12e15]. We also compare this model with other ANN models developed in this work and Machine Learning (ML) classifiers published before to address the same problem. A very good MI-DRAGON QSAR model was obtained, and the subsequent combined QSAR & CN analysis may become of major importance for the prediction of the activity of new compounds against different targets or the discovery of new targets. In this sense we reported an illustrative study that combines both experiment and theory to show how to use this model in practical situations. We reported the prediction, synthesis, and pharmacological assay of oxoisoaporphines with MAO-A inhibitory activity. In Fig. 1 we depict a flowchart with the main steps given in this work to train and validate the ANN classifier (Fig. 2). 2. Materials and methods 2.1. Computational methods 2.1.1. MI-DRAGON technique 2.1.1.1. Parameters for drugs. The DRAGON software 4.0 [11] was utilized here to calculate the parameters of drugs. This software provides near to 1600 descriptors classified as zero- (0D) one- (1D), two- (2D), and three-dimensional (3D) descriptors. It depends on whether they are computed from the chemical formula, substructure list representation, molecular graph or geometrical representation of the molecule, respectively [16,17]. In this work, we calculated only 0D eto e2D descriptors. We use these descriptors because this way the drugs not optimized for use with 3D descriptors. We used specifically the following descriptors: 2D autocorrelations, Burden eigenvalues, topological charge indices, Fig. 2. Snapshot of LOMETS server used to predict 3D structure of peptides. F. Prado-Prado et al. / European Journal of Medicinal Chemistry 46 (2011) 5838e58515840 eigenvalue-based indices, functional group counts, atoms-centered fragments, charge descriptors and molecular properties. 2.1.1.2. Parameter for targets (proteins). In previous works, we have predicted protein function based on 3D parameters of protein structure calculated with software MI [18,19]. In this paper, we used specifically the 3D-electrostatic potential parameters x k . These values were used as inputs to construct the QSAR model together with the structural parameters of drugs. The detailed explanation of the procedure has been published before. Therefore, we provide herein only the more general formula for these potentials and some general explanations [20]: x m¼ x kðRÞ¼ X n j¼1˛R pkðjÞ$ x ðjÞ(1) The average general potentials depend on the absolute probabilities p k (j) and the total potential with which the aminoacid j-th interact with the rest of aminoacids. These are the probabilities with which the amino acids interact with other amino acids placed at a distance equals to k-times the cut-off distance (r ij ¼k$r cut-off ). The method uses an MCM to calculate these probabilities; which also depend on the 3D interactions between all pairs of aminoacids placed at distance r ij in r 3 in the protein structure. However, for the sake of simplicity, a truncation or cut-off function a ij is applied in such a way that a short-term interaction takes place in a first approximation only between neighboring aa ( a ij ¼1ifr ij <r cut-off ). Otherwise, the interaction is banished ( a ij ¼0). The relationship a ij may be visualized in the form of a protein structure complex network. In this network the nodes are the C a atoms of the aminoacids and the edges connect pairs of aminoacids with a ij ¼1. This network can be understood in terms of aminoacideaminoacid protein contact maps [21] in Euclidean 3D space r 3 ¼(x,y,z) coordinates of the C a atoms of aminoacids listed on protein PDB files. In recent works we published different examples of these networks [22,23]. For the purposes of the calculation, all water molecules and metal ions were removed [20]. All calculations werecarried out with our in-house software MI [20]. For calculation the MI software always uses the full matrix, never a sub-matrix, but may run the last summation term either for all amino acids or only for some specific groups called regions (R). These regions are often defined in geometric terms and are referred to as core, inner, middle and surface regions. The protein regions (c correspond to core, i to inner, m to middle, and s to surface regions, respectively) are shown in different figures published in previous works [24]. The diameters of the regions, as a percentage of the longest distance r max with respect to the centre of charge, are 0e25 for region c, 26 to 50 for region i, 51 to 75 for region m, and 76 to 100 for regions. Additionally, we consider the total region (t) that contains all the amino acids in the protein (region diameter 0e100% of r max ). Consequently, we can calculate different x m values for a single protein. Each x m values is referred for all the amino acids contained in different Regions R of the protein (c, i, m, s, or t) and all their neighbors placed at topological distance k withinthis region (kis the namedthe order). In this work, we calculated in total 5 regions 6 (higher order considered) ¼30 x m values for each protein. 2.1.1.3. Statistical analysis. The linear mt-QSAR model was constructed using Linear Discriminant Analysis (LDA). All statistical analyses and data exploration were carried out in STATISTICA 6.0 [25]. In so doing, we used the Forward stepwise method implemented in the LDA module of STATISTICA for the selection of variables. In the present work, the independent data test is used by splitting the data randomly in a training series used for a model construction and a cross-validation (CV) one. Let be d k drugs and x m target molecular descriptors of type kth for different drugs (d), we attempt to develop a simple linear classifier of mt-QSAR type with the general formula: SðDTPÞpred ¼X ndd k¼0 ak$dkþX ntd m¼0 bm$ x mþc0(2) We used LDA to fit this discriminant function. The model deals with the classification of a compound set with or without affinity on different receptors. A dummy variable Affinity Class (AC) was used as input to codify the affinity. This variable indicates either high (AC ¼1) or low (AC ¼0) affinity of the drug by the receptor. S(DPI) pred or DPI affinity predicted score is the output of the model and it is a continuous dimensionless score that sorts compounds from low to high affinity to the target coinciding DPIs with higher values of S(DPI) pred and nDPIs with lowest values. In Equation (2), aandbrepresents the coefficients of the classification function, determined by the LDA module of the STATISTICA 6.0 software package [25]. We used Forward Stepwise algorithm for a variable selection. The statistical significance of the LDA model was determined calculating the p-level (p) of error with Chi-square test. We also inspected the Specificity, Sensitivity, and total Accuracy to determine the quality-of-fit to data in training. The validation of the model was corroborated with external prediction series. 2.1.1.4. ANN analysis. The non-linear mt-QSAR model was constructed using ANN analysis. All models trained were carried out in STATISTICA 6.0 [25]. In so doing, we used a very simple type of ANN called Three Layers Perceptron (MLP-3) to fit this discriminant function. The model deals with the classification of a compound set with or without affinity on different receptors. A dummy variable Affinity Class (AC) was used as input to codify the affinity. This variable indicates either high (AC ¼1) or low (AC ¼0) affinity of the drug by the receptor. S(DTP) pred or DTP affinity predicted score is the output of the model and it is a continuous dimensionless score that sorts compounds from low to high affinity to the target coinciding DTPs with higher values of S(DTP) pred and nDTPs with lowest values. In Equation (2),brepresents the coefficients of the LNN classification function, determined by the ANN module of the STATISTICA 6.0 software package [25]. We used Forward Stepwise algorithm for a variable selection. In addition, we can explore more complicated non-linear ANNs in order to improve the accuracy of the classifier. We processed our data with different ANNs looking for a better model. Four types of ANNs were used, namely, Probabilistic Neural Network (PNN), Radial Basic Function (RBF), Linear Neural Network (LNN), and Four Layer Perceptron (MLP-4) [26,27]. The quality of all the ANNs (linear or non-linear) was determined calculating values of Specificity, Sensitivity, and total Accuracy to determine the quality-of-fit to data in training. The validation of the model was corroborated with external prediction series. We also reported ROC-curve analysis (ROC curve can be used to select an optimum decision) for both training and validation series [26,28]. 2.1.1.5. Data set. The data set was formed by a set of marketed DPIs with known affinity of drugs by targets. This dataset is the same benchmark data used in previous works [2,8e10] in this area and contains all drugs approved by the US FDA. We download this dataset from the public resource called Drug Bank [8,9]{Knox, 2011 #11246; Wishart, 2008 #11249; Wishart, 2006 #11250}. The data set was formed for more than 519 drugs with their respectively 336 targets. Subsequently, we were able to collect above 2337 cases (drug-protein interactions) instead of 519 336 cases. In addition the data set was used to develop ANN models to performance the model. The names or codes for all compounds are depicted in F. Prado-Prado et al. / European Journal of Medicinal Chemistry 46 (2011) 5838e5851 5841 Table 1SM and Table 2SM of the supplementary material, due to space constraints, as well as the references consulted to compile the data in this table. 2.1.1.6. Complex network construction. We construct a DPIs network in order to achieve the drug and protein affinity with a network approach. Generally in this network, one node may represent a drug or a target. On the other hand, the edges represents the DPIs; express relationships between pairs of drugs with their targets [10]. Anyhow, the nodes representing targets may be of at least two types. In almost all cases reported up to date each target is represented only once in the network. In this class of “static”DPIs network the target is depicted by the node corresponding to the X-ray structure of itself. In this work, we build in total two complex networks. First, we constructed the DPIs networks for the observed data and second, DPIs network predicted by the model. The common steps to construct these networks are: 1. First, using the Excel software in a column we introduce all the proteins, the drugs used quotation marks in our database. 2. Then in another column lists all the cases. At the beginning of this column puts the total number of vertices, there are currently two columns of the name of drug and protein and their corresponding number of vertices. 3. At the end of the columns are placed bows in the first column put the number of vertices for the drug and in another column corresponding to the protein. 4. The file was saved as a .txt format file. After we had renamed the .txt file as a .net file we read it with the CentiBin software [29,30]. 5. Using CentiBin we can not only represent the network but also highlight all drugs and targets (nodes) connected by a specific edge or link (DPI). Using this software we can calculate vertex centralities to analyze the relationships between drug targets. 2.2. Illustrative experiments 2.2.1. Synthesis of oxoisoaporphines 2.2.1.1. Synthesis. Synthesis of compounds 1e10 has been previously reported by us [31,32]. The 2,3-dihydrooxoisoaporphines 2 and 3 (Scheme 1) wereobtained starting from condensationproduct of 3,4-dimethoxyphenylethylamine (homoveratrylamine) with phthalaldehydic acid. The intermediate compound was subsequently treated with polyphosphoric acid to give the compounds 2 and 3. Treatment with 10% Pd on charcoal over benzene of 2 and 3 yielded the isoaporphines 5 and 6 39. Compound 2 was catalytically hydrogenated over PtO2 at room temperature at 60e70 psi in AcOH affording, in good yield, 9 in which only ring D is saturated 40, 41. By NaBH4 reduction of 2 carbinol 10 is obtained 37. Finally, treatment of 5 with dust zinc and 37% HCl to give the phenolic compound 7. On the other hand, 1 (Scheme 2) was obtained from N-phenethylphthalimide, which was partially reduced with NaBH 4 in MeOH and cyclized with hydrochloric acid to give 5,6,8,12b-tetrahydro-8isoindolo [1,2-a]isoquinolone and this was oxidized with air in the presence of NaOH-MeOH and dimethyl sulfate to afford 1-(2methoxycarbonylphenyl)-3,4-dihydroisoquinoline, which was directly hydrolyzed with hydrochloric acid to 1-(2-carboxyphenyl)- 3,4-dihydroisoquinoline, which, using fuming sulfuric acid, was finally cyclized affording 1. Treatment of 1 with 10% Pd on charcoal over benzene yielded isoaporphine 4 [33]. The BischlereNapieralski condensation of homoveratrylamine and 2-(benzylbenzoate)chloride afforded a compound characterized as (20-(3,4-dihydro-6,7dimethoxyisoquinolin-10-yl)phenyl)methylbenzoate that, when it was made to react with an AcOH/H 2 SO 4 mixture at 100  C, surprisingly afforded only compound 5-methoxy-6H-dibenzo[de,h] quinolin-6-one compound 8 (Scheme 3) [34,35]; Schemes 1, 2 and 3 are in Fig. 3. 3. Results 3.1. DPIs QSAR predictive models 3.1.1. LDA model Common physicochemical properties like electrostatic potentials have been demonstrated to be useful on protein QSAR [36,37]. We used these properties as input of our model in addition to drug molecular descriptors. The present is the first mt-QSAR model combining DRAGON and MI to predict the probability with which occur DPIs between a drug and a protein. This type of models lie within the frontiers between classic QSAR for drugs and protein QSAR [38]. Some applications for the present model are the prediction of new drugs, new protein receptors or drug targets, and drug binding sites. Detailed information on the compounds, predicted classification, and probability of affinity on different receptors of the drugs used to seek the model appears in Table 1SM of the supplementary material. Based on the algorithms described in materials and methods the best linear model found was the following: SðDTPÞpred ¼0:0004$d10:0032$d20:0026$d3 þ0:0021$d40:0019$d50:0003$d6 0:0026$d7þ0:1284$d80:0253$ x 1þ4:2642 N¼2337 c 2¼2585:53 plevel <0:001 ð3Þ The nomenclature used in the descriptors of the equation is found in Table 1. In this equation, Nis the number of cases, c 2 is the Chi-square and p is the level of error. This model, with 9 variables, classifies correctly 533 out of 588 DPIs (Sensitivity of 77.47%) and 800 out of 810 nDPIs (Specificity of 98.76%). Overall training Accuracy was 88.98%. The validation of the model was carried out by means of external predicting series. The model classifies correctly 252 out of 339 DPIs (74.34%) and 490 out of 510 nDPIs (96.08%) in validation series. Accuracy for validation series (predictability) was 87.37%. These results (Table 2) indicate that we developed an accurate model according to previous reports on the use of LDA in QSAR [39,40]. After observing the above result, we developed another model to incorporate possible interaction effects between drug and protein descriptors, with the aim of improving the model. For this, we used the product between descriptors resulting from the linear mt-QSAR equation. In fact, one interaction effect was entered in the new equation after rerunning forward stepwise analysis. The model equation incorporating interaction effects is as follows: SðDTPÞpred ¼0:0044$d2þ0:0013$d9þ0:0032$d10 þ0:00022$d60:0017$d11 0:0003$d6 þ0:0015$d4$ x 11:643 N¼2337 c 2¼565:4862 plevel <0:001ð4Þ This model, with 9 variables, classifies correctly 276 out of 339 DPIs (Sensitivity of 81.42%) and 430 out of 510 nDPIs (Specificity of 84.31%). Overall training Accuracy was 83.16%. The validation of the model was carried out by means of external predicting series. The model classifies correctly 539 out of 678 DPIs (79.5%) and 700 out of 810 nDPIs (86.42%) in validation series. Accuracy for validation series (predictability) was 83.27%. These results (Table 2) indicate that we can seek a statistically significant and relatively accurate linear model according to previous reports on the use of LDA in QSAR [39,40]. We can also conclude that even when the stepwise F. Prado-Prado et al. / European Journal of Medicinal Chemistry 46 (2011) 5838e58515842 analysis incorporates an interaction effect it does not improve the final model. In conclusion, there is a not-random linear relationship between DPIs and drug þprotein descriptors calculated with DRAGON and MI. We can also conclude that we may need to carry out non-linear techniques to improve the model. 3.1.2. Train and validation of 2D MI-DRAGON ANN model The previous models show good results with a relatively small number of parameters (9 parameters and 7 parameters) and a linear equation for each one. However, as result of the previous section we decided to carry out an ANN analysis to seek a better model using a non-linear method. Four types of ANNs were used, namely, Probabilistic Neural Network (PNN), Radial Basic Function (RBF), Three Layers Perceptron (MLP-3), and Four Layer Perceptron (MLP-4). See, previous works about the use of these ANNs in protein QSAR [2,9]. The Fig. 4 depicts the networks topology for some of the ANN models tested. In general, at least one ANN of every type tested was statically significant. However, one must note that the profiles of each network indicate that many of these are highly non-linear and complicated models. ANN-QSAR model has been demonstrated before; see, for instance, the works of Fernandez and Caballero [41,42].We Fig. 3. Oxoisoaporphines derivatives used in this work. F. Prado-Prado et al. / European Journal of Medicinal Chemistry 46 (2011) 5838e5851 5843 compare different types of networks to obtain a better model. In Table 2 we show the classification matrix of the different networks. The profiles of networks tested were RBF 1:1-1-1:1 with only one variable; LNN 243:243-1:1, which present many variables, and PNN 243:243-7458-2-2:1, which has a very high number of hidden neurons, see Table 2. After that, the simpler but more accurate ANN model found was an MLP (MLP 21:21-31-1:1) with training Accuracy ¼91.06%. This was selected as the best network found because it presents both high accuracy and an adequate number of variables accounting for features relevant for DPIs. This ANN presents 21 inputs variables (18 d k þ3 x m ). This leads to 21 neurons in first or input layer (I), 31 neurons in the second layer or first hidden layer (H1) and only one neuron (DPI prediction) in the output layer (O). We depict the ROC-curve for MLP 21:21-31-1:1 to show how reliable was the network model developed, see Fig. 5. Notably, almost all the models presented had a ROC-curve higher than 0.5. The model presented an area greater than 0.92. From now on we call the ANN MLP 21:21-31-1:1 as the 2D MI-DRAGON predictor. 3.1.3. 2D MI-DRAGON assembly of CNs for DPIs A possible application for this model, which is relevant to drug and target screening, is the construction of multi-protein CNs that incorporates protein affinity profile for drugs or the same CNs for DPIs. In order to recall the capacity of 2D MI-DRAGON to predict new CNs of DPIs we selected the same benchmark database used in previous works [2,9]; which includes US FDA approved drugs with their targets. With these goals in mind, we constructed again and manually curated the above-mentioned CN obtaining a graph with 855 vertices or nodes (drugs and proteins) and m¼1016 DPIs (edges). This CN of DPIs have D¼6.7; average topological distances D ij between all pairs of nodes. The same as before, we constructed a new CN of DPIs but connecting only pairs of nodes with DPIs predicted by 2D MI-DRAGON. In so doing, we obtained a value of D¼6.6 and m¼907 DPIs. In Fig. 6 we illustrated visually both CNs (observed and predicted). In first instance, we can see that both networks may have very similar topology (connectivity patterns structure) as measure in terms of D and m. One way to apply this type of CNs of DPIs for drug screening and drug-target discovery is the calculation of those nodes (drugs or proteins) which are more relevant or important (central) in the graph. For it we can use numerical parameters that quantify the importance of a node in a graph which are called node centralities C t of type t [43]. The identification of these nodes using node centralities may help us to identify the more relevant drugs or Table 2 Comparison of LDA and different ANNs classification models. Model profile Class Train Stat. Validation % DPIs nDPIs Par. % DPIs nDPIs 2D MI-DRAGON DPIs 89.38 303 36 Sn 92.18 625 53 MLP nDPIs 94.12 30 480 Sp 90.12 80 730 21:21-31-1:1 Total 92.23 Ac 91.06 LDA a DPIs 77.47 155 533 Sn 77.47 155 533 9:9e1:1 nDPIs 98.77 800 10 Sp 98.77 800 10 Total 88.99 Ac 88.99 LDA DPIs 84.31 430 80 Sn 86.42 700 110 6:6e1:1 nDPIs 81.42 63 276 Sp 79.50 139 539 Total 83.16 Ac 83.27 PNN DPIs 63.72 216 123 Sn 67.55 458 220 243:243-7458-2-2:1 nDPIs 74.51 130 380 Sp 67.90 260 550 Total 70.20 Ac 67.74 RBF DPIs 74.04 251 88 Sn 71.24 483 195 1:1-1-1:1 nDPIs 72.55 140 370 Sp 80.25 160 650 Total 73.14 Ac 76.14 LNN DPIs 84.37 286 53 Sn 80.09 543 135 243:243-1:1 nDPIs 100.00 0 510 Sp 81.48 150 660 Total 93.76 Ac 80.85 DPIs: Drug-Target Pairs for compounds with high affinity; nDPIs: Drug-Target Pair for compounds with non-affinity; Stat. is statistics, Par. is parameter. Table 1 Detailed list of the symbols and description for all parameters present in the model. Original Descriptor Descriptor name Code ID ATS6v Broto-Moreau autocorrelation of a topological structure - lag 6/weighted by atomic van der Waals volumes d 1 GATS2e Geary autocorrelation - lag 2/weighted by atomic Sanderson electronegativities d2 BELm1 Lowest eigenvalue n. 1 of Burden matrix/weighted by atomic masses d3 BELm4 Lowest eigenvalue n. 4 of Burden matrix/weighted by atomic masses d4 BELm5 Lowest eigenvalue n. 5 of Burden matrix/weighted by atomic masses d5 GGI1 Topological charge index of order 1 d6 GGI9 Topological charge index of order 9 d7 JGI7*10-3 Mean topological charge index of order7/1000 d8 T q ðmÞEntropy of all aminoacids placed in the middle region and all the neighbors at distance k2p1 GATS6p Geary autocorrelation - lag 6/weighted by atomic polarizabilities d9 BELv5 Lowest eigenvalue n. 5 of Burden matrix/weighted by atomic van der Waals volumes d10 SEigv Eigenvalue sum from van der Waals weighted distance matrix d11 d4p1 (lowest eigenvalue n. 4 of Burden matrix/weighted by atomic masses) (Entropy of all aminoacids placed in the middle region and all the neighbors at distance k2) d4p1 Fig. 4. Generic Topology of ANN models trained in this work. F. Prado-Prado et al. / European Journal of Medicinal Chemistry 46 (2011) 5838e58515844 human Guests: theory, applications, legal Protection, taxes, and regulatory issues, Curr. Proteomics 6 (2009) 214e227. [54] H. González-Díaz, A. Duardo-Sanchez, F.M. Ubeira, F. Prado-Prado, L.G. PérezMontoto, R. Concu, G. Podda, B. Shen, Review of MARCH-INSIDE & complex networks prediction of drugs: ADMET, anti-parasite activity, Metabolizing Enzym. Cardiotoxicity Proteome Biomarkers Curr. Drug Metab. 11 (2010) 379e406. [55] H. Gonzalez-Diaz, F. Romaris, A. Duardo-Sanchez, L.G. Perez-Montoto, F. Prado-Prado, G. Patlewicz, F.M. Ubeira, Predicting drugs and proteins in parasite infections with topological indices of complex networks: theoretical backgrounds, applications, and legal issues, Curr. Pharm. Des 16 (2010) 2737e2764. [56] C. Yang, L.G. Valerio Jr., K.B. Arvidson, Computational toxicology approaches at the US Food and drug Administration, Altern. Lab. Anim. 37 (2009) 523e531. F. Prado-Prado et al. / European Journal of Medicinal Chemistry 46 (2011) 5838e5851 5851 3D MI-DRAGON: new model for reconstruction of US FDA drug-target network and theoretic-experimental studies of rasagiline derivatives inhibitors for AChE. Francisco Prado-Prado 1*, Xerardo García-Mera 1, Manuel Escobar1, Nerea Alonso1, Olga Caamaño1,Matilde Yañez 2, and Humberto González-Díaz 3 1 Department of Organic Chemistry, Faculty of Pharmacy, University of Santiago de Compostela (USC), 15782, Spain. 2 Department of Pharmacology, Faculty of Pharmacy, USC, 15782, Spain 3 Department of Microbiology and Parasitology, Faculty of Pharmacy, USC, 15782, Spain Abstract. The Neurodegenerative diseases have been increasing in the last years. Many of the drugs candidates to be used in the treatment of neurodegenerative disease present specific 3D structural features. One important protein in this sense is the acetylcholinesterase (AChE); which is the target of many Alzheimer's dementia drugs. Consequently, the prediction of Drug-Proteins Interactions (DPIs/nDPIs) between new drugs candidates with specific 3D structure and targets it is of the major importance. For it, we can use Quantitative Structure-Activity Relationships (QSAR) models to carry out rational DPIs prediction. Unfortunately, many previous QSAR models developed to predict DPIs take into consideration only 2D structural information and codify the activity against only one target. To solve this problem we can develop one 3D multi-target QSAR (3D mtQSAR) models. In this communication, we introduce the technique 3D MI-DRAGON a new predictor for DPIs based two different well-known software. We use the software MARCH-INSIDE (MI) and DRAGON to calculate 3D structural parameters for drugs and targets respectively. Both classes of 3D parameters were used as input to train Artificial Neuronal Network (ANN) algorithms using as benchmark dataset the complex network (CN) formed by all DPIs between US FDA approved drugs and their targets. The entire dataset was downloaded from Drug Bank. The best 3D mt-QSAR predictor found is one ANN of type Multi-Layer Perceptron (MLP) with profile MLP 37:37-24-1:1. This MLP classifies correctly 274 out of 321 DPIs (Sensitivity = 85.35%) and 1041 out of 1190 nDPIs (Specificity = 87.48%), corresponding to training Accuracy = 87.03%. We validated the model with external predicting series with Sensitivity = 84.16% (542/644 DPIs; Specificity = 87.51% (2039/2330 nDPIs) and Accuracy = 86.78%. The new CNs of DPIs reconstructed from US FDA can be used to explore large DPIs databases in order to discover both new drugs and/or targets. We carried out theoretic-experimental studies to illustrate the practical use of 3D MI-DRAGON. First, we reported the prediction and pharmacological assay of 22 different rasagiline derivatives with possible AChE inhibitory activity. Keywords: Drug-Protein interaction complex networks; Protein Structure Networks; multi-target QSAR; Markov Model; AChE inhibitors Corresponding authors: PRADO PRADO, F. ([email protected]), Faculty of Pharmacy, University of Santiago de Compostela 15782, Spain., Fax: +34-981 594912. *Manuscript Click here to view linked References 1. Introduction Yildirim, et al. [1] have built a complex network (CN) of Drug-Protein Pairs (DPIs) with the form of a bipartite graph composed of all DPIs for all US Food and Drug Administration (US FDA) approved drugs and proteins linked by drug-target binary associations. The resulting CN connects most drugs into a highly interlinked giant component, with strong local clustering of drugs of similar types according to Anatomical Therapeutic Chemical classification. It was motivated due to the strong incentive to develop new methods able of predicting potential drug-target interactions complex networks (CNs) formed by DPIs [2]. For it, we can use Quantitative Structure-Activity Relationships (QSAR) models [5] to carry DPIs prediction. To solve this problem we can develop a 3D multi-target QSAR (3D mt-QSAR) models to predict DPIs [6]. One way to develop this class mt-QSAR is incorporating into the QSAR equation parameters of the structure of the target (protein, DNA, RNA, etc.) in addition to the structural parameters of the drug present in classic QSAR. Some of the more known software we can use to reach this goal are: DRAGON, CODESSA[7], MODES-LAB[8], TOMO-COMD[9], and MARCH-INSIDE (MI)[10]. The software DRAGON is one of the more complete calculating more than 1600 descriptors for drug structure including as zero- (0D) one- (1D), two- (2D), three-dimensional (3D) parameters. Unfortunately several QSAR models are able to predict the activity of drugs against only one target and/or are unable to codify important 3D structural features. Speck-Planche, et al.[3, 4] have developed mt-QSAR for the design of multi-target inhibitors against chemokine receptors. This approach was focused on the construction of a mt-QSAR model for the classification and prediction of inhibitors chemokine receptors. For instance, very recently we have developed in a previous work a QSAR model base on the MARCHINSIDE method to predict a large network of DTPs [11].This model was based on 2D structural parameters for drugs and 1D structural parameters for protein. After that we developed MIND-BEST [12] and NL MIND-BEST [13]. Both predictors are based on 3D structural parameters of proteins calculated with software MI but they used only 2D structural parameters of drugs (calculated also with MI). The accuracy of the MIND-BEST model found was 86.32% and NL MIND-BEST was Accuracy = 90.41%. However both models only use 2D parameters using MI software. After that, to improve and obtain better results we use the software MARCH-INSIDE (MI) to calculate 3D structural parameters for targets and the software DRAGON was used to calculated 2D molecular descriptors all drugs[14]. We introduce the technique 2D MI-DRAGON a new predictor for DPIs based on two different well-known software. As was mentioned in the previous paragraph we can seek a QSAR predictor for DPIs using molecular descriptors of both drug and target. In this work, we introduce for first time 3D MI-DRAGON a new predictor for DPIs based on two different well-known software. We use the software MARCH-INSIDE (MI) to calculate 3D structural parameters for targets and the software DRAGON for 3D parameters of all DPIs present in the Drug Bank (US FDA benchmark dataset) [15-18]. Both classes of parameters were used as input of different Artificial Neuronal Network (ANN) algorithms to seek an accurate non-linear mtQSAR predictor. 3D MI-DRAGON offers a good opportunity for fast-track calculation of all possible DPIs of one drug enabling us to re-construct large drug-target or DPIs Complex Networks (CNs). In this study, we reported the prediction and pharmacological assay of 22 different rasagiline derivatives with AChE inhibitory activity. The present work reports the attempts to calculate within unified DPIs. All this can help to design new inhibitors of AChE. A very good 3D MI-DRAGON QSAR model was obtained, and the subsequent combined QSAR & CN analysis may become of major importance for the prediction of the activity of new compounds against different targets or the discovery of new targets. In this sense we reported an illustrative study that combines both experiment and theory to show how to use this model in practical situations. We reported the prediction and pharmacological assay of rasagiline derivatives with AChE inhibitory activity. In Figure 1 we depict a flowchart with the main steps given in this work to train and validate the ANN classifier. Figure 1 comes about here 2. Materials and Methods 2.1. Computational methods 2.1.1 MOPAC AM1 Optimization geometry method using CS CHEM 3D. Molecular structures of all FDA drugs were generated with CHEM 3D Ultra (version 2005). The energy of each intermediate was then minimized using the semi-empirical MOPAC method with a minimum RMS gradient of 0.100, which specifies the convergence criteria for the gradient of the potential energy surface. The geometry of the molecules was optimized and the values of the quantum chemical descriptors of each compound were calculated using AM1. AM1 theory was used with a closed shell function. The MOPAC AM1 method was selected because it was a semi-empirical quantum chemical method and the computational time was much shorter than that needed by ab initio method. 2.1.2. MI-DRAGON technique 3D Parameters for drugs. The DRAGON software 4.0 [19] was utilized here to calculate the 3D parameters of drugs. It depends on whether they are computed from the chemical formula, substructure list representation, molecular graph or geometrical representation of the molecule, respectively [20, 21]. In this work, we calculated only GETAWAY 3D descriptors. We use these descriptors after optimized for use with 3D descriptors. 3D Parameters of proteins. In previous works we have predicted protein function based on different protein structural parameters derived from a Markov matrix that account for electrostatic interactions between aminoacid pairs in the 3D structure of the protein. One of the classes of parameters used was called the Shannon Entropy Tθk(R) of the markov matrix. These values are used here as inputs to describe information about the structure of the drug target proteins (T) in order to construct the mt-QSAR models for DTPs. The detailed explanation has been published before [22-30] and reviewed in detail more recently [31]. At follows we give the formula for Tθk(R) values and some general explanations:           1log    Rj j k j k k TRpRpR  Where, kpi(R) values are the absolute probabilities with which the effect of the electrostatic interaction propagates from the amino acid ith to other amino acids jth next to it and returns to ith after k-steps. These probabilities refer to: aminoacids considered isolated in the space (k = 0), interaction between aminoacids in direct contact (k = 1) or spatial (k > 1) indirect interactions between amino acids placed at a distance equal to k-times the cut-off distance (rij = k ·rcut-off) in the residue network. Euclidean 3D space r3 = (x, y, z) coordinates of the Cα atoms of amino acids listed in protein PDB files. For calculation, all water molecules and metal ions were removed [32]. All calculations were carried out with our inhouse software MARCH-INSIDE 2.0 [32]. For the calculation, the MARCH-INSIDE software always uses the full matrix, never a sub-matrix, but the last summation term may run either for all amino acids or only for some specific protein regions (R) denoted as: c for core, i for inner, m for middle, and s for surface regions, respectively). Consequently, we can calculate different Tθk(R) for the amino acids contained in the regions (c, i, m, s, or t) and placed at a topological distance k each other within this orbit (k is the order) [22, 23, 33-35]. In this work, we have calculated altogether 5(types of regions) x 6(orders considered) = 30 Tθk(R) indices for each protein. 2.1.3 Statistical analysis. Let be Dθk(G) entropy descriptors molecular that codify information about drug structure and Tθk(R) entropy descriptors that codify information about drug target proteins; we attempt to develop a simple mt-QSAR model in the form of a linear classifier with the general formula:         2 0 5 0 , 5 0 ,cRbGaDTPS k k T kR k k D kG pred     We used Linear Discriminating Analysis (LDA) to fit this discriminant function. The model deals with the classification of a compound set with or without affinity on different receptors. A dummy variable Affinity Class (AC) was used as input to codify the affinity. This variable indicates either high (AC = 1) or low (AC = 0) affinity of the drug by the receptor. S(DTP)pred or DTP affinity predicted score is the output of the model and it is a continuous dimensionless score that sorts compounds from low to high affinity to the target coinciding DTPs with higher values of S(DTP)pred and nDTPs with lowest values. In equation (6), b represents the coefficients of the classification function, determined by the LDA module of the STATISTICA 6.0 software package [36]. We used Forward Stepwise algorithm for a variable selection. The statistical significance of the LDA model was determined calculating the p-level (p) of error with Chi-square test. We also inspected the Specificity, Sensitivity, and total Accuracy to determine the quality-of-fit to data in training. Cases for training set were selected at random out of the cases in full dataset. The remnant cases were used to validate the model. The validation of the model was corroborated with these external prediction series; these cases were never used to train the model. The ration between training/validation set was 2/1 approximately. This procedure to select training and validation sets is largely known and used to train QSAR models [37-43]. 2.1.4 ANN analysis. The non-linear mt-QSAR model was constructed using ANN analysis. All models trained were carried out in STATISTICA 6.0 [36]. In so doing, we used a very simple type of ANN called Three Layers Perceptron (MLP-3) to fit this discriminant function. The model deals with the classification of a compound set with or without affinity on different receptors. A dummy variable Affinity Class (AC) was used as input to codify the affinity. This variable indicates either high (AC = 1) or low (AC = 0) affinity of the drug by the receptor. S(DTP)pred or DTP affinity predicted score is the output of the model and it is a continuous dimensionless score that sorts compounds from low to high affinity to the target coinciding DTPs with higher values of S(DTP)pred and nDTPs with lowest values. In equation (2), b represents the coefficients of the LNN classification function, determined by the ANN module of the STATISTICA 6.0 software package [36]. We used Forward Stepwise algorithm for a variable selection. In addition, we can explore more complicated non-linear ANNs in order to improve the accuracy of the classifier. We processed our data with different ANNs looking for a better model. Four types of ANNs were used, namely, Probabilistic Neural Network (PNN), Radial Basic Function (RBF), Linear Neural Network (LNN), and Four Layer Perceptron (MLP-4)[44, 45]. The quality of all the ANNs (linear or non-linear) was determined calculating values of Specificity, Sensitivity, and total Accuracy to determine the qualityof-fit to data in training. The validation of the model was corroborated with external prediction series. We also reported ROC-curve analysis (ROC curve can be used to select an optimum decision) for both training and validation series [44, 46]. 2.1.5 Data set. The data set was formed by a set of marketed DPIs with known affinity of drugs by targets. This dataset is the same benchmark data used in previous works [1, 5, 12, 13] in this area and contains all drugs approved by the US FDA. We download this dataset from the public resource called Drug Bank [12, 13, 16-18]. The data set was formed for more than 519 drugs with their respectively 336 targets. Subsequently, we were able to collect above 4485 cases (drug-protein interactions) instead of 519 x 336 cases. In addition the data set was used to develop ANN models to performance the model. The names or codes for all compounds are depicted in Table 1SM of the supplementary material, due to space constraints, as well as the references consulted to compile the data in this table. 2.1.6 Complex network construction. We construct a DPIs network in order to achieve the drug and protein affinity with a network approach. Generally in this network, one node may represent a drug or a target. On the other hand, the edges represents the DPIs; express relationships between pairs of drugs with their targets [1]. Anyhow, the nodes representing targets may be of at least two types. In almost all cases reported up to date each target is represented only once in the network. In this class of “static” DPIs network the target is depicted by the node corresponding to the X-ray structure of itself. In this work, we build in total two complex networks. First, we constructed the DPIs networks for the observed data and second, DPIs network predicted by the model. The common steps to construct these networks are: First, using the Excel software in a column we introduce all the proteins, the drugs used quotation marks in our database. Then in another column lists all the cases. At the beginning of this column puts the total number of vertices, there are currently two columns of the name of drug and protein and their corresponding number of vertices. After, at the end of the columns are placed bows in the first column put the number of vertices for the drug and in another column corresponding to the protein. Then, the file was saved as a .txt format file. After we had renamed the .txt file as a .net file we read it with the CentiBin software [47, 48]. Finally, using CentiBin we can not only represent the network but also highlight all drugs and targets (nodes) connected by a specific edge or link (DPI). Using this software we can calculate vertex centralities to analyze the relationships between drug targets. 2.2. Illustrative experiments 2.2.1. Synthesis of Rasagiline derivatives. Synthesis. Synthesis of compounds 1-22 has been previously reported by us [5, 12], see Figure 2. Figure 2 comes about here 2.2.2 Determinations of cholinesterases activities The cholinesterase assay method of Ellman was used to determine the in vitro cholinesterase activity [49]. The activity was measured by increase in absorbance at 412 nm due to the yellow color produced from the reaction of acetylthiocholine iodide with the dithiobisnitrobenzoate (DTNB) ion. Acetylcholinesterase from human erythrocytes, acetylcholinesterase recombinant expressed in HEK 293 cells and butyrylcholinesterase from human serum was obtained from Sigma. 2.2.3 Experimental conditions and kinetics. Enzyme activity was measured using a FLUOstar Optima microplate reader. The assay medium contained phosphate buffer, pH 8.0, 20 mM DTNB, 0.01 U/ml of enzyme and 0.75 µM substrate (acetylthiocholine iodide or butyrylthiocholine iodide). The activity was determined by measuring the increase in absorbance at 412 nm at 1 min intervals for 10 min at 37 °C. In dose-dependent inhibition studies, the substrate was added to the assay medium containing enzyme, buffer, and DTNB with inhibitor after 10 min of incubation time. All experiments were carried out in duplicate and expressed as mean ± SEM. The relative activity is expressed as percentage ratio of enzyme activity in the absence of inhibitor, see Table 1. Table 1 comes about here 3. Results 3.1. DPIs QSAR predictive models 3.1.1 LDA model. Common physicochemical properties like entropy have been demonstrated to be useful on protein QSAR [50, 51]. We used these properties as input of our model in addition to drug molecular descriptors. The present is the first mt-QSAR model combining DRAGON and MI to predict the probability with which occur DPIs between a drug and a protein. This type of models lie within the frontiers between classic QSAR for drugs and protein QSAR [33]. Some applications for the present model are the prediction of new drugs, new protein receptors or drug targets, and drug binding sites. Detailed information on the compounds, predicted classification, and probability of affinity on different receptors of the drugs used to seek the model appears in Table 1SM of the supplementary material. Based on the algorithms described in materials and methods the best linear model found was the following:   001.0919.29884485 )3(48.225.010.011.065.1 62.117.5234.977.1237.3601.11 2 10987 654321    levelpN dddd ddddddDTPSpred  Table 2 comes about here The nomenclature used in the descriptors of the equation is found in Table 2. In this equation, N is the number of cases, χ2 is the Chi-square and p is the level of error. This model, with 10 variables, classifies correctly 256 out of 321 DPIs (Sensitivity of 79.75%) and 1014 out of 1190 nDPIs (Specificity of 85.21%). Overall training Accuracy was 84.05%. The validation of the model was carried out by means of external predicting series. The model classifies correctly 498 out of 644 DPIs (77.33%) and 2000 out of 2330 nDPIs (85.84%) in validation series. Accuracy for validation series (predictability) was 83.99%. These results (Table 3) indicate that we developed an accurate model according to previous reports on the use of LDA in QSAR [52, 53]. Table 3 comes about here 3.3.2 3D MI-DRAGON ANN model. The previous model show good results with a relatively small number of parameters (10 parameters) and a linear equation. However, as result of the previous section we decided to carry out an ANN analysis to seek a better model using a non-linear method. Four types of ANNs were used, namely, Probabilistic Neural Network (PNN), Radial Basic Function (RBF), Three Layers Perceptron (MLP-3), and Four Layer Perceptron (MLP-4). See, previous works about the use of these ANNs in protein QSAR [5, 13]. The Figure 3 depicts the networks topology for some of the ANN models tested. In general, at least one ANN of every type tested was statically significant. However, one must note that the profiles of each network indicate that many of these are highly non-linear and complicated models. Figure 3 comes about here Models using ANN-QSAR has been demonstrated before; see, for instance, the works of Fernandez and Caballero [54, 55]. We compare different types of networks to obtain a better model. In Table 3 we show the classification matrix of the different networks. The profiles of networks tested were RBF 1:1-1-1:1 with only one variable; LNN 227:227-1:1, which present many variables, and PNN 227:227-14797-2-2:1, which has a very high number of hidden neurons, see Table 3. After that, the simpler but more accurate ANN model found was an MLP (MLP 37:37-24-1:1) with training Accuracy = 87.03 %. This was selected as the best network found because it presents both high accuracy and an adequate number of variables accounting for features relevant for DPIs. This ANN presents 37 inputs variables (24 dk + 13 Θm). This leads to 37 neurons in first or input layer (I), 24 neurons in the second layer or first hidden layer (H1) and only one neuron (DPI prediction) in the output layer (O). We depict the ROC-curve for MLP 37:37-24-1:1 to show how reliable was the network model developed, see Figure 4. Notably, the model presented had a ROC curve higher than 0.5. The model presented an area greater than 0.92. From now on we call the ANN MLP 37:37-24-1:1 as the 3D MI DRAGON predictor. Figure 4 comes about here 3.3.2.1 3D MI-DRAGON assembly of CNs for DPIs., The construction of multi-protein CNs that incorporates protein affinity profile for drugs or the same CNs for DPIs is relevant to drug and target screening. And is one application for this model. In order to recall the capacity of 3D MI-DRAGON to predict new CNs of DPIs we selected the same benchmark database used in previous works [5, 13, 14]; which includes US FDA approved drugs with their targets. With these goals in mind, we constructed again and manually curated the above-mentioned CN obtaining a graph with 855 vertices or nodes (drugs and proteins) and m = 1016 DPIs (edges). This CN of DPIs have D = 6.7; average topological distances Dij between all pairs of nodes. The same as before, we constructed a new CN of DPIs but connecting only pairs of nodes with DPIs predicted by 3D MI-DRAGON. In so doing, we obtained a value of D = 7.2 and m = 1256 DPIs. In Figure 5 we illustrated visually both CNs (observed and predicted). Figure 5 comes about here In first instance, we compare this predicted network (3D MI-DRAGON) with 2D MIDRAGON predicted network[14]. We compare to observe the similar or dissimilar topology (connectivity patterns structure) between them. Measuring in terms of TIs such as: number of nodes (n), number of edges (m), Wiener index (W), diameter (D), the Randic connectivity index (Xr), topological distance (Dist), network average values for radiality (R), node degree (δ), eccentricity (E). In Table 4, we observe all the TIs are similar excepting n, m and w. That means both CNs have a high similarity between them. These results are very interesting, because our 3D MI-DRAGON model present similar results to the 2D MI-DRAGON model, which results have been published successfully before. Table 4 comes about here To see how reliable and valid is our model. Not only compared to TIs to observe similarity between both predicted networks, but we study the centrality analysis of given networks too. This type of drug screening and drug target discovery is the calculation of those nodes (drugs or proteins) which are more relevant or important (central) in the graph. For it we can use numerical parameters that quantify the importance of a node in a graph which are called node centralities Ct of type t [56]. These nodes identification using node centralities may help us to identify the more relevant drugs or proteins in analogy to similar procedures developed for PINs; networks of Protein-Protein Interactions (PPIs) [57]. In Table 5 we show the predicted results of both node degree centrality (Cδ) and closeness centrality (Cclo) for proteins and drugs present in the database and compare with the predicted results of 2DMI-DRAGON model. The parameter Cδ measures the local importance of a node by counting the number of nodes directly attached to him [57]. Conversely, Cclo measures the global importance of a node in a CN by taking in consideration the inverse of the sum of Dij (Cclo = 1/ΣDij) [58]. Consequently, the higher Cδ the higher is the local importance of the node but the higher Cclo the lower is the global importance of the node. For instance, the protein 1HA2 is one important protein both locally and globally in this CNs with lower Cclo > 4 and a Cδ = 26. It means that this protein is both locally and globally important because it is the target of many drugs (high Cδ). This result is similar to obtained by 2D MIDRAGON model. Another interesting result was simvastatin. Simvastatin is a hypolipidemic drug used to control elevated cholesterol, or hypercholesterolemia. It is a member of the statin class of pharmaceuticals. The primary use of simvastatin is for the treatment of dyslipidemia and the prevention of cardiovascular disease [59, 60]. Depending on our aims the more important nodes in pharmacological terms not necessarily have to be the more central in the graph (those with higher Cδ and lower Cclo), see Table 5. We show in this example, our model predicts efficiently. We found that the 3D MI-DRAGON model shows very similar results to the previous model, which has been published with excellent results. In Table 2SM of supplementary material we show all node degrees and closeness results. Table 5 comes about here 3.2. Theoretic-Experimental Study using 3D MI-DRAGON predictor Finally, we illustrated in one theoretic-experimental study the practical use of 3D MIDRAGON. We reported the prediction, synthesis, and pharmacological assay of 20 different rasagiline derivatives with AChE inhibitory activity. 3.2.1 3D MI-DRAGON prediction of rasagiline derivatives vs. AChE. In this in silico experiment we used 3D MI-DRAGON to predict the interaction of the rasagiline derivatives with respect to AChE. For it, we downloaded the 3D structure of AChE protein with PDB ID 1EEA and calculated their structural parameters with MI. We also generated the SMILE codes for these compounds and we use MOPAC AM1 Optimization geometry method for these compounds for calculated their 3D structural parameters with DRAGON. After that, we predicted their propensity to undergo DPIs with AChE using as inputs for the 3D MI-DRAGON predictor the structural parameters of both the drugs and the protein. In Table 6 we confront the results obtained using this model and the outcomes of the pharmacological assay. No compounds are selective inhibitors of AChE, which is why we used as control galantamine for AChE was. We consider the observed class for active compounds OC = 1 if compound IC50 < 10 μM this cutoff is in the similar range than other used in previous works [61, 62]. As we can see in this table all the compounds rasagiline derivatives present some activity, But none of these compounds have inhibitory activity in the pharmacological assays. All of our compounds in the pharmacological assay (OC = 0) were inactive. 3D MI-DRAGON predicted as inactive all compounds, excepting 3. The model classified correctly 19 of 22 compounds tested (86.36%). In this test, our model was compared with pharmacological testing of 22 compounds synthesized by us. And we can observe the effectiveness of our model with experimental data. Also, we note that the model predicts all compounds tested as inactive, this is important because the model allows to discriminate between active and inactive compounds. However, some compounds were not tested by pharmacological assay, that compounds were predicted as inactive using 3D MI-DRAGON model. We discarded pharmaceutical assays of these compounds; because we consider our model reliable. This kind of model can be used to saves efforts and money to perform the pharmacological tests. This is a good example of how reliable is the MI DRAGON 3D model. Table 6 comes about here 3.2.2 3D MI-DRAGON complex network of rasagiline derivatives vs. US FDA proteins. An additional use of 3D MI-DRAGON was to carry out the “in silico” or virtual screening of the new compounds with respect to all other targets previously approved by US FDA [14, 63]. It may help to found new targets for these drugs or discard possible toxicological effects depending on the other targets predicted and/or discarded for these compounds. This type of experiment is of the major importance due to the cost in terms of animal sacrifice, time, materials and human resources of the experimental assay of all compounds against all these targets, see recent reviews by Duardo-Sanchez et al. [64-67]. In fact, over a decade, the US FDA has been engaged in the applied research, development, and evaluation of computational toxicology methods used to support the safety evaluation of a diverse set of regulated products. The basis for evaluating computational toxicology methods is multi-factorial, including the potential for increased efficiency, reduction in the numbers of animals used, lower costs, and the need to explore emerging technologies that support the goals of the US FDA's Critical Path Initiative (e.g. to make decision support information available early in the drug review process)[68]. In this experiment, we downloaded the 3D structure of all proteins that are targets of US FDA approved drugs. Next, we calculated the structural parameters of all these proteins with MI. We also generated the SMILE codes for these compounds and and we use MOPAC AM1 Optimization geometry method for these compounds for calculated their 3D structural parameters with DRAGON. After that, we predicted their propensity to undergo DPIs with all US FDA proteins using as inputs for the 3D MI-DRAGON predictor the structural parameters of both the drugs and proteins. We predicted all proteins in FDA dataset vs. the 22 rasagiline derivatives. We found that most of 22 derivatives were predicted as non-active (low DPIs scores) against most proteins in the FDA database. Consequently, 3D MI-DRAGON predicts a high selectivity of rasagiline derivatives as AChE inhibitors. We can reach this goal because the model predicts these compounds as non-active with respect to most proteins that are targets of FDA drugs. Using these results, we constructed a DP-CN for rasagiline derivatives and the FDA dataset (see Figure 6). As a result we obtained a CN with 87 nodes (FDA drugs, proteins, or rasagiline derivatives) and 166 DP (edges, DTPs). As In this network we can see that protein 1EEA (AChE) is predicted to interacts with compound 3, this protein is an AChE target [69]. These results are good because they agree with the experimental results presented in this paper where the compound 3 show low AChE activity. The use of such complex networks can help us find and predict new drugs-protein interactions, and therefore find new drugs with improved biological activity and fewer side effects, especially in neural disease. Figure 6 comes about here 4. Conclusions Figure Legends: Figure 1. Flowchart of all steps given in this work to develop the new model Figure 2. Rasagiline derivatives used in this work Figure 3. Generic Topology of ANN models trained in this work Figure 4. ROC Curve for 3D MI-DRAGONGON predictor (red = train series, blue = validation series Figure 5. Observed vs. Predicted drug-target complex networks Figure 6. Complex network of rasagiline derivatives vs. US FDA proteins Tables: Table 1. Inhibitory activity of different rasagiline derivatives . Compounds hAChE (IC50  M) hAChE (IC50  M) 1 >100 µM 12 >100 µM 2 >100 µM 13 No tested 3 >100 µM 14 No tested 4 No tested 15 No tested 5 ** 16 No tested 6 No tested 17 ** 7 ** 18 ** 8 No tested 19 ** 9 ** 20 ** 10 >100 µM 21 ** 11 >100 µM 22 ** Galantamine 1.43 ± 0.03a Eserine 151.40 ± 5.63 nM Tacrine 130,90 ± 6,83 nM Each IC50 value is the mean ± S.E.M. from five experiments. Table 2. Detailed list of the symbols and description for all parameters present in the model. Original Descriptor Descriptor name Code ID H7v H autocorrelation of lag 7 / weighted by atomic van der Waals volumes d1 HATS5v leverage-weighted autocorrelation of lag 5 / weighted by atomic van der Waals volumes d2 HATS4e leverage-weighted autocorrelation of lag 4 / weighted by atomic Sanderson electronegativities d3 HATS6e leverage-weighted autocorrelation of lag 6 / weighted by atomic Sanderson electronegativities d4 R5e+ R maximal autocorrelation of lag 5 / weighted by atomic Sanderson electronegativities d5   core T 4  Entropy of all aminoacids placed in the core region and all the neighbors at distance k ≤ 4 d6   core T 5  Entropy of all aminoacids placed in the core region and all the neighbors at distance k ≤ 5 d7   inner T 5  Entropy of all aminoacids placed in the inner region and all the neighbors at distance k ≤ 5 d8   middle T 2  Entropy of all aminoacids placed in the middle region and all the neighbors at distance k ≤ 2 p1   surface T 0  Entropy of all aminoacids placed in the surface region and all the neighbors at distance k ≤ 0 d9 Table 3. Comparison of LDA and different ANNs classification models. Model Train Stat. Validation profile Class % DPIs nDPIs Par. % DPIs nDPIs MI DRAGON 3D DPIs 85.36 274 47 Sn 84.16 542 102 MLP nDPIs 87.48 149 1041 Sp 87.51 291 2039 37:37-24-1:1 Total 87.03 Ac 86.79 LDAa DPIs 79.75 256 65 Sn 77.33 498 146 10:10-1:1 nDPIs 85.21 176 1014 Sp 85.84 330 2000 Total 84.05 Ac 83.99 PNN DPIs 0 0 644 Sn 0 0 321 227:227-14797-2-2:1 nDPIs 100 0 2346 Sp 100 0 1174 Total 78.46 Ac 78.53 RBF DPIs 47.05 303 341 Sn 52.65 169 152 1:1-1-1:1 nDPIs 56.01 1032 1314 Sp 54.86 530 644 Total 54.08 Ac 54.38 LNN DPIs 53.73 346 298 Sn 45.79 147 174 227:227-1:1 nDPIs 32.05 1594 752 Sp 31.52 804 370 Total 36.72 Ac 34.58 DPIs: Drug-Target Pairs for compounds with high affinity; nDPIs: Drug-Target Pair for compounds with nonaffinity; Stat. is statistics, Par. is parameter Table 4 Comparison 3D MI-DRAGON versus 2D MI-DRAGON. 2D MI-DRAGON Value TIs Value 3D MI-DRAGON 706 n 59 907 m 631 1826812 W 2057954 18 D 19 255.09 Xr 266.39 2.49 δ 2.44 6.7 Dist 7.2 0.078 E 0.083 12.43 R 11.39 aThe TIs used are: number of nodes (n), number of edges (m), Wiener index (W), diameter (D), the Randic connectivity index (Xr), topological distance (Dist), network average values for radiality (R), node degree (δ), eccentricity (E). Table 5. Results of node degree (Cδ) and closeness centrality (Cclo) for 20 proteins and drugs. Drug/PDB Cδ 2D-MI-DRAGON Cδ 3D MI-DRAGON Drug/PDB Cclo 2D-MI-DRAGON Cclo 3D MI-DRAGON 1HA2 44 26 1HA2 4.80 3.54 1BNA 36 40 Simvastatin 4.22 4.11 NADH 35 33 Gliclazide 4.17 3.48 1R5K 27 29 Saquinavir 4.16 3.46 Simvastatin 18 17 1BNA 4.15 3.46 1EMI 16 21 Cefalotin 4.13 2.98 1CZM 14 17 Atorvastatin 4.09 3.44 1NHZ 14 14 1A8M 4.09 3.75 1MO8 14 14 1XF0 4.07 3.09 1UZF 13 13 Estrone 4.06 2.61 1SQN 13 13 Ketoprofen 4.02 2.86 1T9N 13 12 Testosterone 4.02 2.42 1BYW 11 15 1TZI 4.02 3.77 1VRU 11 10 1KED 4.00 3.16 1E3G 11 10 Captopril 3.99 3.49 Atorvastatin 11 10 Liothyronine 3.97 3.21 1ZNC 10 11 Diflunisal 3.96 3.20 1ODW 10 10 Halothane 3.95 3.11 Pyridoxal Phosphate 9 9 Digitoxin 3.94 3.45 1HWL 9 9 Pyridoxine 3.94 3.08 Table 6. Prediction of rasagiline derivatives with 3D MI-DRAGON predictor DRUG OC PC Score Structure DRUG OC PC Score Structure 1 0 0 0.95 12 0 0 0.63 2 0 0 0.95 13 0 0 1.00 3 0 0 0.86 14 0 0 0.87 4 0 0 0.95 15 0 0 1.00 5 0 1 0.52 16 0 0 0.87 6 0 0 0.95 17 0 0 1.00 7 0 1 0.52 18 0 0 1.00 23 8 0 0 0.88 19 0 0 1.00 9 0 0 0.89 20 0 0 1.00 10 0 0 0.75 21 0 0 0.97 11 0 1 0.27 22 0 0 0.97 OC = Observed class; PC = Predicted class Figure(s) Click here to download high resolution image Figure(s) Click here to download high resolution image technological, and social networks (Bornholdt and Schuster, 2003; Boccaletti et al., 2006;Dehmer and Emmert-Streib, 2009). A network is a set of items, usually called nodes, with connections between them, which are called links or edges (Newman, 2003). The nodes can be atoms, molecules, proteins, nucleic acids, drugs, cells, organisms, parasites, people, words, laws, computers, or any other part of a real system. The edges or links are relationships between the nodes such as chemical bonds, physical interactions, metabolic pathways, pharmacological action, law recurrence, or social ties. There are many different experimental and/or theoretical methods to assign node–node links depending on the type of network we want to construct. Unfortunately, many of these methods are expensive in terms of time or resources. In addition, different methods to link nodes in the same type of network are not totally accurate in such a way that they do not always coincide. For instance, Modha and Singh, intheirwork‘Networkarchitecture of the long-distance pathways in the macaque brain’ (Modha and Singh, 2010) studied the information contained in the ‘Collation of Connectivity data on the Macaque brain’ (CoCoMac) neuroinformatic database in order to construct the most comprehensive long-distance network of the Macaque brain. This database contains 410 anatomical tracing studies, 10,681 connectivity relations and 16,712 mapping relations, and after collation of all connections, a final network of 383 brain regions and 6602 longdistance brain connections that travel through the brain’s white matter were obtained. However, to construct this network, the authors had to solve problems related with the multiplicity of brain maps, divergent nomenclature, boundary uncertainty, different resolutions depending on the work studied. In this context, the development of fast and cheap computational methods in order to collate connectivity information becomes a goal of major importance. One possible solution to this problem is the use of Quantitative Structure–Activity/Property Relationships (QSAR/QSPR) models, which have been traditionally studied in the field of chemoinformatics and are used to predict the biological activity of drugs (QSAR) or physicochemical properties of organic compounds (QSPR) using as input structural parameters of the system under study (Puzyn et al., 2010). In the case of global studies (properties of full system) these parameters are Topological Indices (TIs) derived from the graphical representation of the system (molecule, etc.). On the other hand, we can use node centralities or local TIs of a sub-graph if we want to predict a local property of part of the system (local chemical reactivity, biotransformation of a toxicophore group in a drug, etc.). Currently, the use of QSPR-like models in which the inputs are graph parameters is not limited to the study of molecules and has been extended to other complex systems (Gonza ´lez-Dı ´az and Munteanu, 2010). Specifically, Shannon entropy is one of the most useful parameters used as input in QSAR/QSPR studies to quantify structural information of molecular graphs (Dehmer et al., 2009). In all the above-mentioned cases, Shannon entropy parameters can be used to quantify structural information locally (nodes, edges, paths, clusters, etc.) and/or globally (full graph). In fact, we have used Markov Chain (MC) to calculate Shannon entropies locally or globally within a graph considering all possible branches at different topological distances. The information is quantified in terms of y k (j) values, which are called the Markov– Shannon entropy node centralities of order kth for all jth states (nodes) of a MC associated to the system. This MC is expressed by a Markov or Stochastic matrix ( P 1 ) and represented by a graph of the studied system. The elements of P 1 are the probabilities 1 p ij with which the ith and jth nodes connect each other (there is a physical or functional tie, link, or relationship) within a graph. Using Chapman– Kolmogorov equations it is straightforward to realize the way to calculate y k (j) values for all nodes in a graph. We can use these values directly or sum some of them to obtain total or local entropies (see Section 2). Our group has introduced the software called MARCHINSIDE (Markovian Chemicals In Silico Design), which has become a very useful tool for QSAR/QSPR studies (Gonzalez-Diaz et al., 2010). This software can calculate 1D (sequence), 2D (connectivity in the plane) and 3D (connectivity in the space) MC parameters, including y k (j) values, for many molecular systems. MARCH-INSIDE is able to characterize small molecules (drugs, metabolites, organic compounds), biopolymers (gene sequence, proteins sequence or 3D structure, and RNA secondary structure) and artificial polymers but can perform a limited manage of other complex networks. It happens because MARCH-INSIDE can read, transform into Markov matrix, represent as graph, and calculate entropies for molecular formats (.mol or SMILE .txt files for drugs, .pdb for proteins, or .ct files for RNAs) but it is unable to upload formats of Complex Networks (.mat, .net, .dat, .gml, etc.). Inthiswork,weuseforthefirsttimeQSPR-likemodelsableto assess the quality of the connectivity of new complex networks assembled with information obtained from many sources not totally accurate. The idea is to seek a QSPR-like model that use as input the Fig. 1. General workflow used in this work. P. Riera-Ferna ´ndez et al. / Journal of Theoretical Biology 293 (2012) 174–188 175 y k (j) values for all possible pairs of nodes in a network to decide which pairs of nodes link each other and which ones do not. This class of model will allow us to computationally re-evaluate all the links in any complex network in such a way that we do not have to rely upon experimentation to confirm the existence or not of a link between all pairs of links. Using this model, we should experimentally confirm only those connections predicted by the model with low link score and/or simply remove them from the network depending on the cost/benefit ratio. As a consequence of this aim, we have re-programmed the MARCH-INSIDE application creating a new software able to manage complex networks. The new program is called MI-NODES (MARCH-INSIDE NOde DEScriptors) and is compatible with other software like Pajek or CentiBin, since it is able to read .mat, .net and .dat formats. A very interesting feature of MI-NODESisthatitcanprocessmultiplenetworkswithinafileand calculate both MC global TIs and/or node centralities for all these networks. It is also able to export them in a single file in networkby-network and or node vs. node output formats. In order to illustrate the use of the new method we have carried out 7 experiments. In each experiment we report for the first time new QSPR models, which are useful to re-evaluate connectivity quality of different types of networks. Although very different systems were studied, the same workflow was used (see Fig. 1). In the first experiment we studied the full metabolic pathway networks of four different organisms (bacteria, yeast, nematode, plant). In the second experiment, we studied different biological networks of parasite–host interactions (PHIs). The third experiment consisted of carrying out a study regarding connectivity quality in the CoCoMac cerebral cortex co-activation network (Modha and Singh, 2010). In the fourth experiment we studied a macroscopic landscape parasitism-spreading network for cattle fasciolosis in NW Spain. In the next experiment we illustrate the application of the method to a complex network for all the historic record of 1940–2004 of the entire Financial Law legislation system (legal–social network) in Spain. All these mentioned networks present one class of nodes or are bipartite. We also studied a 5th-partite network (with five classes of nodes) representing different relationships in the world trade of active & intelligent packaging for food industry. This additional network represents the relationships between companies, countries, trade mark products, products uses, and food types. Finally, we carried out an experiment to seek the first QSPR-like model for the US FDA Drug–Target network. This last experiment has methodological particularities because it involves different systems with graduated structural levels. We used as input y k values of drug molecular graph, protein structural networks, and drug–target network. With the advent of the age of complex systems science this study opens a door to a relatively less studied but very important field: the assessment of the connectivity quality in new complex networks. According a recent comprehensive review (Chou, 2011), to develop a useful predictor for a statistical system, the following things often need to be considered: (i) benchmark dataset construction or selection, (ii) formulation of statistical samples, (iii) operating algorithm (or engine), (iv) anticipated accuracy, and (v) web–server establishment. Below, let us elaborate how to deal with these procedures one by one. 2. Materials and methods 2.1. Datasets used In the present study we have selected 7 types of networks, taking into account the data availability, the size of the networks studied (since we are interested in large complex networks) and the level of organization (from molecules to social sciences). Some of them have been used in previous studies and other are presented for the first time in this work. 2.1.1. Metabolic pathway networks Metabolic network data was downloaded directly from Barabasi’s group web (http://www.nd.edu/~networks/resources. htm) as gzipped ASCII file. In this file each number represents a substrate in the metabolic network of corresponding organism. Data-format is: From -To (directed link). The information studied was previously obtained by Jeong et al. from the ‘intermediate metabolism and bioenergetics’ portions of the WIT database and used in order to try to understand the large-scale organization of metabolic networks (Jeong et al., 2000). According to the authors, biochemical reactions described within a WIT database are composed of substrates and enzymes connected by directed links. For each reaction, educts and products were considered as nodes connected to the temporary educt–educt complexes and associated enzymes. Bidirectional reactions were considered separately. For a given organism with Nsubstrates, Eenzymes and Rintermediate complexes the full stoichiometric interactions were compiled into an (NþEþR)X(NþEþR) matrix, generated separately for each of the different organisms. 2.1.2. Parasite–Host complex networks In order to construct the studied networks we have used two sources of information. The first of them is the Interaction Web Database (IWDB) (http://www.nceas.ucsb.edu/interactionweb/ index.html), which contains datasets on species interactions from several communities in different parts of the world. In particular we have used the host–parasite dataset, composed of data belonging to studies about parasites (nematodes, acanthocephalans, cestodes, trematodes, monogeneans, leeches, copepods and branchiurans) and their hosts (fish) from 7 Canadian freshwater systems (Arai and Mudry, 1983;Arthur et al., 1976;Bangham, 1955;Chinniah and Threlfall, 1978;Dechtiar, 1972;Leong and Holmes, 1981). The second source of information is the Global Mammal Parasite Database (GMPD) (http://www.mammalpara sites.org/), a compilation of records of parasites (helminths, protozoa, viruses, bacteria, arthropods and fungi) and their hosts (wild mammals) that have been documented in the published scientific literature (Nunn and Altizer, 2005). In this work we have used the information about ungulates (Order Artiodactyla and Perissodactyla), carnivores and primates. Based on the data obtained from the two databases, we constructed four bipartite networks (Parasite–Fish, Parasite–Ungulates, Parasite–Carnivores and Parasite–Primates) in which the first set of nodes is composed by parasites and the second by hosts, linked if the parasite interacts with the host. 2.1.3. Cerebral Cortex co-activation network The version of the CoCoMac network used in this work consists of 383 hierarchically organized regions spanning cortex, thalamus, and basal ganglia; models the presence of 6602 directed long-distance connections(isthreetimeslargerthan any previously derived brain network) and contains sub-networks corresponding to classic corticocortical, corticosubcortical, and subcortico-subcortical fiber systems (Modha and Singh 2010), see Fig. 2. 2.1.4. Complex network for fasciolosis spreading in NW Spain The dataset reported by Mezo et al. (2008) in a previous work was used by our group to construct a network of farm-to-farm spreading of fasciolosis in cattle for Galicia (NW Spain) in other work (Gonza ´lez-Dı ´az et al., 2010). In this work each farm was considered as a node of the network associated to a Boolean or connectivity matrix Cwith elements C ij (links). As this is a P. Riera-Ferna ´ndez et al. / Journal of Theoretical Biology 293 (2012) 174–188176 symmetric condition, the existence of a connection ij implies the existence of the inverse connection ji. Loops, connections from j-th to the same j-th farm (representing self-infection of animals inside the same farm) were allowed. We place an arc (directed edge) connecting the i-th farm with the j-th farm if they meet the condition given in the Microsoft Excel command (see next equations) that is used to truncate the farm-to-farm distance function (see also Fig. 2). The connectivity of the network C depends on the input parameters: spatial coordinates (x i ,y i ) of the farm (f i ), the altitude of the place (h i ), and the strength of drug treatment (Tr i ¼0, 1, 2, 3) used for pre-existent fasciolosis in this farm (Tr i ¼0 indicates that the disease was not detected). Consequently the matrix Cquantifies the propensity C ij ¼1of the disease to spread between farms immediately after treatment. On the other hand, matrix Lincludes two criteria: the preexistence of a high propensity for disease spreading C ij ¼1 and the experimental confirmation of a high Risk Ratio (RR ij ) of Prevalence After Treatment (PAT j ) of the disease. See the definition of these networks in mathematical terms using Excel functions: L ij ¼ifðANDðC ij ¼1,RR ij 41Þð1Þ RR ij ¼ðPAT i þ1Þ=ðPAT j þ1Þð2Þ C ij ¼if ðORðd ij 4d cutoff n 1 mSumðd ij Þ,d ij ¼0Þ,0,1Þð3Þ d ij ¼0:5 n ðh i þh j Þ n Tr i n Tr j n SQRððx i x j Þ^ 2þðy i y j Þ^ 2Þð4Þ 2.1.5. Law co-recurrence network of the Spanish Financial-Legal system The studied network is built establishing connections between two laws or legal norms (nodes) if the time-lag is less than 1 for the same type of laws. Consequently, law–law links represent the co-recurrence of the Spanish Financial System along time to different norms depending on socio-economical conditions. The Cutoff function for the Spanish financial law recurrence network associated to the matrix Lwith elements L ij is the following (see Fig. 3): L ij ¼if (($D6–G$3) 4$C$3, 0, if ($B6¼G$1, 1, 0)). In this function $B6 and G$1 are Excel references to the column and row containing one-letter codes used to identify the type of financial law approved (node classes). The absolute Excel reference $C$3 points to the cell containing the time-lag cutoff value t off .In addition, $D6 and G$3 are references to the column and row containing the values of the variable time (yy) equal to the year of approbation for a given law or norm. In Fig. 3 we illustrate the Excel sheet for the assembly of this Legal–Social network. 2.1.6. Active & Intelligent Packaging for food industry Here we studied the world trade complex network for active & intelligent packaging in food industry (year 2011). This network interconnects product trade items of five classes (five node classes). The classes of nodes are: Product (PR), Company (CO), Country (CU), Food Type (FT), and product use identified as Packaging Type (PT). The network is 5th-partite in such a way that nodes of a given class are connected with nodes of other classes but nodes of the same class are never connected to each other. This network has been constructed, manually curated, and studied in a previous work. In this previous work the network was assembled using data obtained from several public resources previously compiled and reviewed in another paper (Pereira de Abreu et al., in press). The network created after these two works and used here contains a total of: 222 different products with registered trade mark (222 nodes of class PR), 60 different companies (CO), 15 countries (CU), 33 Food types (FTs) and 29 product types (PT). It makes a total of 359 nodes interconnected by 3868 links. The information encoded by a link depends on the classes of the two nodes interconnected. For instance, if one node belongs to the class CO and the other to the class PR, this link indicates that this company produces and commercializes this product. 2.1.7. US FDA Drug–Target network In this case, the data used to built the drug–target network were obtained from the DrugBank database (http://www.drugbank.ca/), a bioinformatics and cheminformatics resource that combines detailed drug (i.e. chemical, pharmacological and pharmaceutical) Fig. 2. Top: innermost core for the undirected version of CoCoMac network. The innermost core is a central sub-network that is far more tightly integrated than the overall network. Bottom: Geographical maps of Galicia (NW Spain) showing the location of the 275 sampled farms. A—Observed data for network C: the status of infection (empty circles: F. hepatica free and filled circles: F. hepatica infected) and the treatment administered on each farm are shown (blue: none; red: an anthelmintic effective against fluke mature stages and green: a fasciolicide effective against immature and mature stages). B—Observed data for network L: Distribution of farms according to the presence of F. hepatica infection (gray: uninfected; cyan: infected with a within-herd prevalence o25 % and pink: infected with a within-herd prevalence Z25%). P. Riera-Ferna ´ndez et al. / Journal of Theoretical Biology 293 (2012) 174–188 177 data with comprehensive drug target (i.e. sequence, structure, and pathway) information (Wishart et al., 2006,2008;Wishart, 2010; Knox et al., 2011). In particular, 532 drugs approved by the Federal Food and Drug Administration (FDA) and 315 targets were studied. 2.2. Computational methods 2.2.1. Markov–Shannon entropy centralities for nodes In information theory, entropy is a measure of the uncertainty associated with a random variable. The term by itself in this context usually refers to the Shannon entropy, which quantifies, in the sense of an expected value, the information contained in a message, usually in units such as bits. Equivalently, the Shannon entropy is a measure of the average information content one is missing when one does not know the value of the random variable. The concept was introduced by Claude E. Shannon in his 1948 paper ‘‘A Mathematical Theory of Communication’’ (Shannon 1948). In the present work, we construct the classical Markov matrix ( 1 P ) for each network as follows. First, we downloaded from public resources the connectivity matrix Lor obtain the data about the links between the nodes to assemble L(nby n matrix, where nis the number of vertices). Next, the Markov matrix P is built. It contains the vertices probability (p ij ) based on L. The probability matrix is raised to the power k, resulting ( 1 P ) k , and multiplied by the vector of the initial probabilities ( 0 p j ). The resulting vectors contain the absolute probabilities to reach the nodes moving throughout a walk of length kfrom node n i ( k p j ) for each kand are the base for the entropy centrality ( y k ) calculation: k P¼ 0 Pð 1 P Þ k ¼½ k p 1 þ k p 2 þþ k p j ð5Þ y k ¼X k p j log k p j ð6Þ 2.2.2. MI-NODES software for calculation of Markov–Shannon entropies MI-NODES (MARCH-INSIDE NOde DEScriptors) is a GUI Python/wxPython application used for the calculation of a new class centralities/topological indices of nodes, sub-networks, or full networks. Actually, it should be considered as the Fig. 3. Excel calculation sheet to obtain the matrix Lof Spain Law Financial system (A) and MI-NODEs graphical interface (B). P. Riera-Ferna ´ndez et al. / Journal of Theoretical Biology 293 (2012) 174–188178 generalization of the software MARCH-INSIDE to manage any kind of complex networks (this program was originally designed to study drugs, proteins and nucleic acid structures). MI-NODES calculates new types of node Centralities k C c (j) based on Markov normalized node probabilities without removing each node previously to perform calculations. It also calculates Markov generalizations of different topological indices k TI c (G) of class cand power kfor the graph G. The tool is both Pajek and CentiBin compatible because it reads networks in the following formats: .net, .dat and .mat. We depict MI-NODES interface and an example of one of the steps in the assembly of a .mat file in Fig. 3. MI-NODES can calculate the following types of Markov node centralities and topological indices of order k-th: Shannon Entropy, Spectral Moments, Harary numbers, Wiener indices, Gutman topological indices, Schultz topological indices, Broto indices, Balaban indices, Kier–Hall connectivity indices, Randic ´ connectivity indices, Galvez indices and Leverage indices. 2.2.3. Dataset used to construct the model The first step to obtain the dataset for each model was to calculate the y k values for the linked nodes of the networks of each type from K¼0 to 5 using MI-NODES. These are the positive cases. In the second step, the same was done but with the negative cases (nodes that are no linked in the observed network). These negative cases were chosen randomly by MI-NODES (we did not take into account 100% of them because their number is much higher than the number of positive cases and this difference can have negative effects on the statistical analysis). Finally, 75% of the data (chosen randomly) was used for training the model and the remaining 25% for cross-validation. 2.2.4. Linear Discriminant Analysis (LDA) models Once the values of the Markov–Shannon entropies were obtained, we carried out a Linear Discriminant Analysis (LDA) by means of the STATISTICA software (StatSoft.Inc. et al., 2002). LDA is possibly the most common technique used in QSPR/QSAR studies with TIs of molecular graphs, protein and RNA structure networks, and bio-molecular complex networks. Let S(L ij ) be the output variable of a model used to score the quality of the connection between two nodes i-th and j-th (L ij ¼1). We can use LDA to seek a linear equation with coefficients a i ,a j ,a ij and a 0 . These are the coefficients of the TIs used as input (in this case local node centralities) in the QSAR/QSPR model and the independent term. The terms a ik and a jk refer to all nodes that lie within the k-th neighborhood (placed at least at topological distance d¼k)ofi-th or j-th single nodes, respectively. The term a ijk refers to the differences between the neighborhoods of a pair of nodes, which may be connected or not. We can use different statistical parameters to evaluate the statistical significance and validate the goodness-of-fit of LDA equation: n¼number of cases, w 2 ¼Chi-square, p¼the error level, as well as the Accuracy, Specificity, and Sensitivity of both train and external validation series (Hill and Lewicki, 2006). Generally, we write the linear LDA equation with the parameters mentioned above in the following form (example using the entropy values as TIs), see also Fig. 1: SðL ij Þ¼ X 5 k¼0 a ik y k ðiÞþ X 5 k¼0 a jk y k ðjÞþ X 5 k¼0 a ijk ½ y k ðiÞ y k ðjÞþa 0 ð7Þ In statistical prediction, the following three cross-validation methods are often used to examine a predictor for its effectiveness in practical application: independent dataset test, subsampling test, and jackknife test (Chou and Zhang, 1995). However, of the three test methods, the jackknife test is deemed the most objective (Chou and Shen, 2008). The reasons are as follows. (i) For the independent dataset test, although all the proteins used to test the predictor are outside the training dataset used to train it so as to exclude the ‘memory’ effect or bias, the way of how to select the independent proteins to test the predictor could be quite arbitrary unless the number of independent proteins is sufficiently large. This kind of arbitrariness might result in completely different conclusions. For instance, a predictor achieving a higher success rate than the other predictor for a given independent testing dataset might fail to keep so when tested by another independent testing dataset (Chou and Zhang, 1995). (ii) For the subsampling test, the concrete procedure usually used in literatures is the 5-fold, 7-fold or 10-fold cross-validation. The problem with this kind of subsampling test is that the number of possible selections in dividing a benchmark dataset is an astronomical figure even for a very simple dataset, as elucidated by Chou and Shen (2008) and demonstrated by Eqs. 28–30 in Chou (2011). Therefore, in any actual subsampling cross-validation tests, only an extremely small fraction of the possible selections are taken into account. Since different selections will always lead to different results even for a same benchmark dataset and a same predictor, the subsampling test cannot avoid the arbitrariness either. A test method unable to yield a unique outcome cannot be deemed as a good one. (iii) In the jackknife test, all the proteins in the benchmark dataset will be singled out one-by-one and tested by the predictor trained by the remaining protein samples. During the process of jackknifing, both the training dataset and testing dataset are actually open, and each protein sample will be in turn moved between the two. The jackknife test can exclude the ‘memory’ effect. Also, the arbitrariness problem as mentioned above for the independent dataset test and subsampling test can be avoided because the outcome obtained by the jackknife crossvalidation is always unique for a given benchmark dataset. Accordingly, the jackknife test has been increasingly and widely used by those investigators who have strong math background to examine the quality of various predictors (see, e.g., (Chen et al., 2009;Chou and Shen, 2010;Chou et al., 2011;Gu et al., 2010;Lin et al., 2011;Mohabatkar 2010;Wang et al., 2011;Xiao et al., 2011;Zakeri et al., 2011;Zeng et al., 2009;Zhang et al., 2011)). However, to reduce the computational time, we adopted the independent testing dataset cross-validation in this study as done by many investigators with support vector machine (SVM) as the prediction engine. 2.2.5. Graphical representation and description of the networks Using graphical/diagrammatic approaches to study complicated systems can provide an intuitive picture or useful insights to help in analyzing complicated mechanisms in these systems, as demonstrated by many previous studies on a series of important biological topics, such as enzyme-catalyzed reactions (Andraos, 2008;Chou, 1989;Chou and Forsen, 1980;Zhou and Deng, 1984), protein folding kinetics and folding rates (Chou, 1990), inhibition of HIV-1 reverse transcriptase (Althaus et al., 1993a;Althaus et al., 1993), inhibition kinetics of processive nucleic acid polymerases and nucleases (Chou et al., 1994), drug metabolism systems (Chou 2010), analysis of DNA sequence (Xie and Mo, 2011), and protein sequence evolution (Wu et al., 2010). Recently, the wenxiang diagrams (Chou et al., 1997) were also used to investigate protein–protein interactions (Zhou, 2011a,2011b). In order to characterize the studied systems, we have represented graphically and calculated some parameters of both the observed and reconstructed networks. The information to construct them was obtained from the connectivity matrices in the case of the observed networks and from the output of the LDA models in the case of the reconstructed networks. CentiBin (http://centibin.ipk-gatersleben.de/index.php) (v.1.4.3) P. Riera-Ferna ´ndez et al. / Journal of Theoretical Biology 293 (2012) 174–188 179 (Junker et al., 2006) was used to prepare the networks using the command: Utils -‘prepare for centralities-undirected’. This command performs the following steps: removal of existing loops, removal of parallel edges, reduction to the giant component and transformation into an undirected network. Once the networks were prepared, it was possible to calculate the graph diameter, Wiener index and average distance. The density (taking into account if the graph is unipartite or bipartite), average degree and Randic ´index were also calculated, but using Pajek (v.1.26) (http://vlado.fmf.uni-lj.si/pub/networks/pajek/)(Batagelj and Mrvar 1998;De Nooy et al., 2005). This program was used to represent graphically the prepared versions of the observed and reconstructed networks. 3. Results and discussion As can be seen in the next sub-sections, the models developed using structural information encoded by y k values offers good results in terms of overall accuracy (ranging from 72.3% to 97.2%) for all the studied networks. However, there are important differences between the best and the worst result (24.9%). This fact could suggest that the capacity to encode information by y k is strongly influenced by the structure of the network (similar values for y k in positive and negative cases could result in a low discrimination power) or that the linear model is less suitable for some types of networks. Anyway, more studies are needed to measure the influence of different factors on the performance of the models based on y k and other TIs. The use of each model developed in this section is restricted to networks with the same features than the networks used to develop the models. For example, the metabolic pathway model (constructed taking into account four model organisms) could be used to evaluate the equivalent metabolic networks of other organisms. An interesting possibility for future research would be to study the same networks with other TIs to see if the results improve (for example models with an accuracy of 70%) or to seek models in which various TIs are combined (in this case each type of TI could encode a different type of structural information). 3.1. Model 1: metabolic pathway networks Study of metabolic networks is of great interest in biology because many applications are directly built on the use of cellular metabolism. Biotechnologists modify the cells and use them as cellular factories to produce antibiotics, industrial enzymes, antibodies, etc. In biomedicine, it is possible to cure metabolic diseases through a better understanding of the metabolic mechanisms, and to control infections by making use of the metabolic differences between human beings and pathogens (Rosa da Silva et al., 2008). For example, the network topologybased approach has been used to uncover shared mechanisms in the study of disease comorbidity (Lee et al., 2008). In a cell or microorganism, metabolic pathways are seamlessly integrated through a complex network of cellular constituents and reactions. However, despite the key role of these networks in sustaining cellular functions, their large-scale structure is not well known. Jeong et al. (2000) showed that, despite significant variation in their individual constituents and pathways, these metabolic networks have the same topological scaling properties and show striking similarities to the inherent organization of complex nonbiological systems. An interesting question to answer is: given a regulatory pathway system consisting of a set of proteins, can we predict which pathway class it belongs to? In this sense, Huang et al. (2011) developed a computational method for the classification and analysis of regulatory pathways using graph property, biochemical and physicochemical property, and functional property. In another work carried out by the same author (Huang et al., 2010), a computational method was developed for the analysis and prediction of the metabolic stability of proteins based on their sequential features, subcellular locations and interaction networks. Many pathways are not totally confirmed experimentally but have been computationally deduced using protein or gene alignment techniques. The idea follows more or less the following scheme: similar proteome -similar enzymes -similar metabolome. On the other hand, the experimental determination of the full metabolome including each metabolite and metabolite bio-transformation pathways is not always an easy task. All this aspects determine the necessity of alignmentfree techniques to assess network connectivity quality in existing models of metabolic pathway networks. Here we developed a model to re-evaluate connectivity using as inputs the y k values for nodes in already-known metabolic networks. For this analysis we have used metabolic networks of four model organisms belonging to different domains of the tree of life. These organisms are: Escherichia coli (EC), Saccharomyces cerevisiae (SC), Caenorhabditis elegans (CE), and Oryza sativa (OS). E. coli is a gram negative bacterium that is commonly found in the lower intestine of warm-blooded organisms. Most EC strains are harmless, but some serotypes can cause serious food poisoning in humans. From the point of view of research, it is one of the best studied bacteria, especially in the areas of genetics, biochemistry and metabolism, and many researchers have based their studies on existing metabolic networks of this prokaryotic model (Baldazzi et al., 2010; Costa et al., 2010; Gerlee et al., 2009;Fowler et al., 2009; Konig et al., 2006;Imielinski et al., 2006;Shi et al., 1999;Lin et al., 2005;Ghim et al., 2005;Schmid et al., 2004;Light and Kraulis, 2004;Burgard and Maranas, 2001;Edwards and Palsson, 2000). S. cerevisiae (a species of yeast) is a fungus with industrial importance used as a model for understanding and engineering eukaryotic cell function. In fact, it was the first eukaryotic genome that was fully sequenced, annotated, and made available publicly (Goffeau, 1997). C. elegans is a free-living nematode that has become a popular model for genetic and molecular research, since it is easy to maintain and has a very fast life-cycle (Burglin et al., 1998). It was the first multi-cellular organism to have its genome completely sequenced (Consortium TCeS et al., 1998). In the field of parasitology, comparison between CE and other parasitic nematodes is an interesting method for studying the function and regulation of some parasite genes (Bird and Opperman, 1998). Other interesting feature is that CE is sensitive to the majority of anti-helmintic drugs that are used against parasitic worm infections of humans and livestock. This has provided the opportunity to use molecular genetic techniques in the worm for mode of action studies (Holden-Dye and Walker, 2007). Finally, O. sativa, commonly known as rice, is a plant of the family Poaceae with great economic importance for the human being. Over recent years, it has gained importance as a model organism for genetic and molecular studies. This is due to its relatively small genome of 420 Mb, whose full sequence was released as early as 2002. The many tools and experimental approaches now available for rice have made it the most widely studied model for cereals (Muller and Grossniklaus, 2010). The best model found was SðL ij Þ¼159:16 y 3 ðe i Þ120:70 y 1 ðp j Þ95:42½ y 5 ðe i Þ y 5 ðp j Þ0:26 n¼74,999 w 2 ¼26,093 po0:001 ð8Þ In this equation, S(L ij ) is a real-valued output variable that scores the propensity of the ith input or educt (e i ) (reactant or substrate) to undergo a metabolic transformation into the product (p j ). The parameter y 1 quantifies the information related to P. Riera-Ferna ´ndez et al. / Journal of Theoretical Biology 293 (2012) 174–188180 the position of the input or reactant metabolite and their direct neighbors (k¼1) in the metabolic network. The parameter y 5 quantifies the information related to middle-long range subsequent metabolic transformations of all the neighbors of the product metabolite (k¼5) in the metabolic network. As we can see in the previous equation the w 2 ¼26,093 statistic corresponds toap-levelo0.001, which indicates a significant discrimination between known metabolic reactions and those metabolite transformations which are not experimentally observed. The model presents very good values of Accuracy, Sensitivity, and Specificity for the recognition of links both in training and external validation series see Table 1. In Table 2 we carry out a graphical and numerical comparison of the giant components of the metabolic networks already known vs. those obtained after re-evaluating link quality with our model. 3.2. Model 2: Parasite–Host networks Due to the importance for the human and animal health and therefore for the economy, much attention has been focused on parasite–host interactions (PHIs). The study of these interactions can help us to understand the role of phylogenetic and ecological factors on the parasite–host specificity (Desdevises et al., 2002; Detwiler and Janovy, 2008;Poulin et al., 2011) and to know how parasites affect the ecosystem functioning (Hatcher et al., 2006; Price et al., 1986;Anderson and May, 1979). In this sense, network theory is a useful tool for analyzing this type of interactions (Poulin, 2010). However, the high experimental difficulty inherent to the in situ accurate determination of PHIs makes the possibility of curate PHIs networks using a computational model very interesting. In this work, we used y k to seek a QSPR-like model able to score the quality of PHIs in known networks. The best model found was SðL ij Þ¼82:62½ y 5 ðp i Þ y 5 ðh j Þ5:52 n¼49,218 w 2 ¼21,728 po0:001 ð9Þ In this equation, S(L ij ) is a real-valued output variable that scores the propensity of the ith parasite specie (p i ) to infect a given host specie (h j ). The Chi-square statistic ( w 2 ) presents a low p-level o0.001, which indicates a significant discrimination between well-established host–parasite relationships and not confirmed parasitism. The model presents good values of Accuracy, Sensitivity, and Specificity for the recognition of parasite– host relationships (links) both in training and external validation series (see Table 1). Consequently, with this simple linear model we could re-evaluate connectivity quality in the already-known PHIs networks in a fast and non-expensive way (without to experimentally resample all PHIs in the correspondent ecological niche). The giant components of the observed and reconstructed networks are described numerically and graphically in Table 3. 3.3. Model 3: Cerebral Cortex co-activation network Connectivity is the key to understanding distributed and cooperative brain functions. Detailed and comprehensive data on large-scale connectivity between primate brain areas have been collated systematically from published reports of experimental tracing studies (Kotter, 2004). Databasing the brain’s anatomical connectivity as delivered by tracing studies is of particular importance as these data characterize fundamental Table 1 Training and Cross-validation results for all models developed in this work. QSPR Training series Model Cross-validation series model NL L % Parameters % NL L 1Metabolic pathway networks 46,029 18,490 NL 71.3 Specificity 71.5 NL 15,384 6,123 2,295 8,185 L78.1 Sensitivity 77.8 L775 2,719 Total 72.3 Accuracy 72.4 Total 2Parasite-host networks 42,576 2,052 NL 95.4 Specificity 95.6 NL 14,144 652 1,275 3,315 L72.2 Sensitivity 71.3 L436 1,085 Total 93.2 Accuracy 93.3 Total 3Cerebral Cortex co-activation network 31,886 2,698 NL 92.2 Specificity 92.5 NL 10,637 867 1,425 3,527 L71.2 Sensitivity 70.4 L488 1,162 Total 89.6 Accuracy 89.7 Total 4NW Spain Fasciolosis Landscape-Spreading network 18,153 149 NL 99.2 Specificity 99.1 NL 6,068 58 405 964 L70.4 Sensitivity 74.2 L116 334 Total 97.2 Accuracy 97.4 Total 5Legal–social network of the Spanish financial law system 15,564 3,401 NL 82.1 Specificity 81.6 NL 5,172 1,169 0 14,986 L100.0 Sensitivity 100.0 L0 5,014 Total 90.0 Accuracy 89.7 Total 6World Trade Intelligent-active food packaging network 27,387 1,623 NL 94.4 Specificity 94.4 NL 9,128 542 692 2,209 L76.1 Sensitivity 77.5 L218 749 Total 92.7 Accuracy 92.9 Total 7US FDA Drug-target network 3,206 981 NL 76.6 Specificity 77.8 NL 1,079 308 189 532 L73.8 Sensitivity 70.4 L71 169 Total 76.2 Accuracy 76.7 Total Rows: Observed classifications; Columns: Predicted classifications; L: Linked; NL: Not linked. P. Riera-Ferna ´ndez et al. / Journal of Theoretical Biology 293 (2012) 174–188 181 structural constraints of the complex and poorly understood functional interactions between the components of real neural systems. The eventual impact and success of connectivity databases, however, will require the resolution of several methodological problems that currently limit their use. These problems comprise four main points: (i) objective representation of coordinate-free, parcellation-based data, (ii) assessment of the reliability and precision of individual data, especially in the presence of contradictory reports, (iii) data mining and integration of large sets of partially redundant and contradictory data, and (iv) automatic and reproducible transformation of data between incongruent brain maps (Stephan et al., 2001). In order to address points ii and iv, we have developed a specific model for the ‘collation of connectivity data on the macaque brain’ (CoCoMac) database (http://www.cocomac.org). The best model found was SðL ij Þ¼70:56 y 1 ðiÞþ74:51 y 5 ðjÞ1:75 n¼39,536 w 2 ¼22,249 po0:001 ð10Þ In this equation, S(L ij ) is a real-valued output variable that scores the propensity of the ith cerebral cortex region to undergo co-activation with the jth region in the CoCoMac network. Table 2 Comparison of observed vs. re-constructed metabolic pathway networks (giant components). Observed network Network descriptors a Reconstructed network Caenorhabditis elegans 1,173 n 1,037 2,842 m 2,214 4.85 Ad 4.27 0.00413 den 0.00412 479.26 R 371.21 6,304,312 W 4,868,792 14 D 14 4.59 AD 4.53 Escherichia coli 2,268 n 1,991 5,620 m 4,393 4.96 Ad 4.41 0.00219 den 0.00222 871.81 R 646.79 22,914,068 W 17,450,524 12 D 12 4.46 AD 4.40 Oryza sativa 658 n 566 1,498 m 1,185 4.55 Ad 4.19 0.00693 den 0.00741 282.24 R 221.01 2,061,882 W 1,494,724 14 D 12 4.77 AD 4.67 Sacharomices cerevisae 1511 n 1295 3807 m 2928 5.04 Ad 4.52 0.00334 den 0.00349 600.49 R 433.95 10,332,580 W 7,455,256 14 D 14 4.53 AD 4.45 b The size of each node is proportional to its normalized degree. a Network descriptors: Total number of connected nodes (n), number of edges (m), average degree (Ad), density (den), Randic ´connectivity index (R), Wiener index (W), diameter (D) and average distance (AD). P. Riera-Ferna ´ndez et al. / Journal of Theoretical Biology 293 (2012) 174–188182 The parameter y 1 (i) quantifies the information related to the position of the ith region and their direct neighbors (k¼1) in the network. The parameter y 5 (j) quantifies the information related to middle-long range co-activation of different brain areas (k¼5) in the cerebral cortex. As in the previous equation the w 2 ¼22,249 statistics corresponds to a p-level o0.001, which indicates a significant discrimination between co-activated regions and not co-activated ones. The model presents very good values of Accuracy, Sensitivity, and Specificity (see Table 1). The giant components of the observed and reconstructed networks are described numerically and graphically in Table 4. 3.4. Model 4: Fasciolosis spreading network (NW Spain) Fasciolosis is a parasitic infection caused by Fasciola hepatica (liver fluke) that has become an important cause of lost productivity in livestock worldwide. Considered a secondary zoonotic disease until the mid-1990s, human fasciolosis is at present Table 3 Comparison of observed vs. re-constructed parasite-host interaction networks (giant components). Observed network Network descriptors a Re-constructed network Parasites–fish 298 n 271 239 np 233 59 nh 38 912 m 788 6.12 Ad 5.82 0.0647 den 0.0890 95.77 R 80.91 298,626 W 246,480 7D6 3.37 AD 3.37 Parasites–Ungulates 793 n 675 701 np 645 92 nh 30 1863 m 1393 4.70 Ad 4.13 0.0289 den 0.0720 191.02 R 125.79 2,534,140 W 1,763,848 10 D 6 4.03 AD 3.88 Parasites–carnivores 619 n 537 537 np 505 82 nh 32 1343 m 1074 4.34 Ad 4.00 0.0305 den 0.0665 159.16 R 115.30 1,587,048 W 1,179,176 9D8 4.15 AD 4.10 Parasites–primates 913 n 612 757 np 582 156 nh 30 1,993 m 1145 4.37 Ad 3.74 0.0169 den 0.0656 252.66 R 118.58 3,609,008 W 1,450,756 10 D 8 4.33 AD 3.88 a Network descriptors: Total number of connected nodes (n), number of connected parasites (np), number of connected hosts (nh), number of edges (m), average degree (Ad), density (den), Randic ´connectivity index (R), Wiener index (W), diameter (D) and average distance (AD). b The size of each node is proportional to its normalized degree. c Outer nodes: Parasites; Inner nodes: Hosts. P. Riera-Ferna ´ndez et al. / Journal of Theoretical Biology 293 (2012) 174–188 183 emerging or re-emerging in many countries, including increases of prevalence and intensity and geographical expansion. In fact, research in recent years has justified the inclusion of fasciolosis in the list of important human parasitic diseases. At present, fasciolosis is the vector-borne disease presenting the widest latitudinal, longitudinal and altitudinal distribution known. In addition, it presents a range of epidemiological characteristics related to a wide diversity of environments (Mas-Coma, 2005). In this sense, the study of geographical spreading of fasciolosis becomes a subject of great interest. In fact, in a recent work we have constructed a network to study the landscape spreading of fasciolosis in Galicia (NW Spain) (Gonza ´lez-Dı ´az et al., 2010). However, we do not have quantitative criteria on the quality of the network connectivity, and re-sampling of all data to reevaluate this connectivity in a field study is a hard and expensive task in terms of time and resources. This situation has prompted us to seek a model in order to assess the quality of the network previously assembled. The best QSPR model found was SðL ij Þ¼20:23 y 1 ðf i Þþ165:13 y 4 ðf j Þ0:82 n¼19,671 w 2 ¼16,058 po0:001 ð11Þ The entropy values y k (f i ) and y k (f j ) used in this equation quantify information about the connectivity patterns between farms in the network C. As can be seen in the equations described in Section 2, the connectivity of Cdepends on the spatial coordinates (x i ,y i ) of the farm (f i ), the altitude of the place (h i ), and the anti-parasite drug treatment (Tr j ) used to prevent Fasciolosis in this farm. Consequently the matrix Cquantifies the a priori propensity C ij ¼1 of this disease to spread between farms immediately after treatment depending on geographical conditions. On the other hand, matrix Lincludes both criteria: (i) the preexistence of a high propensity for disease spreading C ij ¼1 and (ii) the experimental confirmation L ij ¼1 of a high Risk Table 4 Comparison of observed vs. re-constructed networks (giant components). Observed network Network descriptors a Reconstructed network Fasciolosis Landscape-Spreading Network for NW Spain 259 n 252 1,798 m 1,278 13.88 Ad 10.14 0.0538 den 0.0404 114.47 R 91.50 273,596 W 260,884 11 D 10 4.09 AD 4.12 Cerebral Cortex co-activation Network (CoCoMac) 360 n 351 5208 m 3,720 28.93 Ad 21.20 0.0806 den 0.0606 145.98 R 115.90 292,892 W 284,238 5D4 2.27 AD 2.31 Financial Law Network 223 n 223 16,611 m 16,611 148.98 Ad 148.98 0.671 den 0.671 111.12 R 111.12 65,790 W 65,790 2D2 1.33 AD 1.33 US FDA Drug–Target Network 638 n 448 404 nd 355 234 nt 93 802 m 527 2.51 Ad 2.35 0.00848 den 0.0160 213.33 R 133.41 2,702,728 W 1,449,898 17 D 18 6.65 AD 7.24 a Network descriptors: Total number of connected nodes (n), number of connected drugs (nd), number of connected targets (nt), number of edges (m), average degree (Ad), density (den), Randic ´connectivity index (R), Wiener index (W), diameter (D) and average distance (AD). b The size of each node is proportional to its normalized degree. c In the Drug–Target network the outer nodes represent drugs and the inner nodes targets. P. Riera-Ferna ´ndez et al. / Journal of Theoretical Biology 293 (2012) 174–188184