draft for ANVIL impl

src/mol-model-props/computed/topology.ts

@@ -6,12 +6,9 @@
  */
 
 import { ParamDefinition as PD } from '../../mol-util/param-definition';
-import { Structure, Unit } from '../../mol-model/structure';
+import { Structure } from '../../mol-model/structure';
 import { CustomStructureProperty } from '../common/custom-structure-property';
 import { CustomProperty } from '../common/custom-property';
-import { QuerySymbolRuntime } from '../../mol-script/runtime/query/compiler';
-import { CustomPropSymbol } from '../../mol-script/language/symbol';
-import Type from '../../mol-script/language/type';
 import { CustomPropertyDescriptor } from '../../mol-model/custom-property';
 import { ANVILParams, Topology } from './topology/ANVIL';
 
@@ -21,22 +18,6 @@ export const TopologyParams = {
 export type TopologyParams = typeof TopologyParams
 export type TopologyProps = PD.Values<TopologyParams>
 
-export const TopologySymbols = {
-    // TODO is this delegation needed?
-    isMembrane: QuerySymbolRuntime.Dynamic(CustomPropSymbol('computed', 'topology.is-membrane', Type.Bool),
-        ctx => {
-            if (!Unit.isAtomic(ctx.element.unit)) return false;
-            return !TopologyProvider.get(ctx.element.structure).value;
-        }
-    ),
-    isNotMembrane: QuerySymbolRuntime.Dynamic(CustomPropSymbol('computed', 'topology.is-not-membrane', Type.Bool),
-        ctx => {
-            if (!Unit.isAtomic(ctx.element.unit)) return false;
-            return TopologyProvider.get(ctx.element.structure).value;
-        }
-    ),
-};
-
 export type TopologyValue = Map<number, Topology>
 
 export const TopologyProvider: CustomStructureProperty.Provider<TopologyParams, Topology> = CustomStructureProperty.createProvider({
@@ -50,7 +31,7 @@ export const TopologyProvider: CustomStructureProperty.Provider<TopologyParams,
     getParams: (data: Structure) => TopologyParams,
     isApplicable: (data: Structure) => true, 
     // TODO needs ASA to be computed (or 'resolved' before trying computing topology) - how to achieve?
-    // TODO potentially, this should behave like secondary structure info where data can be either parsed or computed
+    // TODO potentially, this could behave like secondary structure info where data can be either parsed or computed
     obtain: async (ctx: CustomProperty.Context, data: Structure, props: Partial<TopologyProps>) => {
         const p = { ...PD.getDefaultValues(TopologyParams), ...props };
         return { value: await Topology.compute(data, p).runInContext(ctx.runtime) };

src/mol-model-props/computed/topology/ANVIL.ts

@@ -29,8 +29,7 @@ export type ANVILProps = PD.Values<ANVILParams>
 export { Topology };
 
 interface Topology {
-    readonly serialResidueIndex: ArrayLike<number>
-    readonly exposure: ArrayLike<Topology>
+    readonly membrane: Vec3[]
 }
 
 namespace Topology {
@@ -53,12 +52,13 @@ namespace Topology {
     const vec = Vec3();
     export async function calculate(runtime: RuntimeContext, structure: Structure, params: ANVILProps): Promise<Topology> {
         const { label_atom_id, x, y, z } = StructureProperties.atom;
-        const elementCount = structure.elementCount;
+        const { label_comp_id } = StructureProperties.residue;
+        const elementCount = structure.polymerResidueCount;
         centroidHelper.reset();
         l.structure = structure;
 
-        const offsets = new Int32Array(elementCount);
-        const exposed: Array<Boolean> = new Array(elementCount);
+        let offsets = new Int32Array(elementCount);
+        let exposed: boolean[] = new Array<boolean>(elementCount);
 
         // ensure ASA
         const accessibleSurfaceArea = structure && AccessibleSurfaceAreaProvider.get(structure);
@@ -86,7 +86,6 @@ namespace Topology {
 
                 // while iterating use first pass to compute centroid
                 Vec3.set(vec, x(l), y(l), z(l));
-                // console.log(vec);
                 centroidHelper.includeStep(vec);
 
                 // keep track of offsets and exposed state to reuse
@@ -97,49 +96,265 @@ namespace Topology {
             }
         }
 
-        centroidHelper.finishedIncludeStep();
-        console.log(centroidHelper.center);
+        // omit potentially empty tail
+        offsets = offsets.slice(0, m);
+        exposed = exposed.slice(0, m);
 
-        for (let k = 0; k < m; k++) {
+        // calculate centroid and extent
+        centroidHelper.finishedIncludeStep();
+        const centroid = centroidHelper.center;
+        for (let k = 0, kl = offsets.length; k < kl; k++) {
             setLocation(l, structure, offsets[k]);
             Vec3.set(vec, x(l), y(l), z(l));
-            // console.log(vec);
             centroidHelper.radiusStep(vec);
         }
+        const extent = 1.2 * Math.sqrt(centroidHelper.radiusSq);
+
+        const initialHphobHphil = HphobHphil.filtered(offsets, exposed, structure, label_comp_id, () => true);
+        const initialMembrane = findMembrane(generateSpherePoints(params.numberOfSpherePoints, centroid, extent), centroid, params, initialHphobHphil, offsets, exposed, structure, label_comp_id);
+        const alternativeMembrane = findMembrane(findProximateAxes(initialMembrane, params.numberOfSpherePoints, centroid, extent), centroid, params, initialHphobHphil, offsets, exposed, structure, label_comp_id);
+
+        const membrane = initialMembrane.qmax! > alternativeMembrane.qmax! ? initialMembrane : alternativeMembrane;
 
-        console.log(1.2 * Math.sqrt(centroidHelper.radiusSq));
-        
         return {
-            exposure: new Array(),
-            serialResidueIndex: new Array()
+            membrane: createMembraneLayers(membrane, extent, params)
         };
     }
 
-    function setLocation(l: StructureElement.Location, structure: Structure, serialIndex: number) {
-        l.structure = structure;
-        l.unit = structure.units[structure.serialMapping.unitIndices[serialIndex]];
-        l.element = structure.serialMapping.elementIndices[serialIndex];
-        return l;
+    function createMembraneLayers(membrane: Membrane, extent: number, params: ANVILProps): Vec3[] {
+        const out: Vec3[] = [];
+        const radius = extent * extent;
+        const normalVector = membrane.normalVector;
+
+        createMembraneLayer(out, membrane.planePoint1, normalVector, params.membranePointDensity, radius);
+        createMembraneLayer(out, membrane.planePoint2, normalVector, params.membranePointDensity, radius);
+        
+        return out;
+    }
+
+    function createMembraneLayer(out: Vec3[], point: Vec3, normalVector: Vec3, density: number, radius: number) {
+        const d = -Vec3.dot(normalVector, point);
+        for (let i = -1000, il = 1000; i < il; i += density) {
+            for (let j = -1000, jl = 1000; j < jl; j += density) {
+                const rep = Vec3.create(i, j, -(d + i * normalVector[0] + j * normalVector[1]) / normalVector[2]);
+                if (Vec3.squaredDistance(rep, point) < radius) {
+                    out.push(rep);
+                }
+            }
+        }
+    }
+
+    interface Membrane {
+        planePoint1: Vec3,
+        planePoint2: Vec3,
+        stats: HphobHphil,
+        normalVector: Vec3,
+        point: Vec3 | undefined,
+        qmax: number | undefined,
+        membraneAtoms: Vec3[] | undefined
+    }
+
+    namespace Membrane {
+        const out = Vec3();
+        export function initial(c1: Vec3, c2: Vec3, stats: HphobHphil): Membrane {
+            Vec3.sub(out, c1, c2);
+            return {
+                planePoint1: c1,
+                planePoint2: c2,
+                stats: stats,
+                normalVector: out,
+                point: void 0,
+                qmax: void 0,
+                membraneAtoms: void 0
+            }
+        }
+
+        export function scored(c1: Vec3, c2: Vec3, stats: HphobHphil, spherePoint: Vec3, qmax: number): Membrane {
+            Vec3.sub(out, c1, c2);
+            return {
+                planePoint1: c1,
+                planePoint2: c2,
+                stats: stats,
+                normalVector: out,
+                point: spherePoint,
+                qmax: qmax,
+                membraneAtoms: []
+            }
+        }
+    }
+
+    function findMembrane(spherePoints: Vec3[], centroid: Vec3, params: ANVILProps, initialStats: HphobHphil, offsets: ArrayLike<number>, exposed: ArrayLike<boolean>, structure: Structure, label_comp_id: StructureElement.Property<string>): Membrane {
+        // best performing membrane
+        let membrane: Membrane | undefined = void 0;
+        // score of the best performing membrane
+        let qmax = 0;
+
+        // construct slices of thickness 1.0 along the axis connecting the centroid and the spherePoint
+        const diam = Vec3();
+        const normedDiam = Vec3();
+        for (let i = 0, il = spherePoints.length; i < il; i++) {
+            const spherePoint = spherePoints[i];
+            Vec3.sub(diam, centroid, spherePoint);
+            Vec3.scale(diam, diam, 2);
+            const diamNorm = Vec3.magnitude(diam);
+            Vec3.scale(normedDiam, diam, 1 / diamNorm);
+            const qvartemp = [];
+
+            const c1 = Vec3();
+            const c2 = Vec3();
+            for (let i = 0, il = diamNorm - params.stepSize; i < il; i += params.stepSize) {
+                const norm = Vec3.magnitude(diam);
+                Vec3.scaleAndAdd(c1, spherePoint, diam, i / norm);
+                Vec3.scaleAndAdd(c2, spherePoint, diam, (i + params.stepSize) / norm);
+
+                // evaluate how well this membrane slice embeddeds the peculiar residues
+                const stats = HphobHphil.filtered(offsets, exposed, structure, label_comp_id, (testPoint: Vec3) => isInMembranePlane(testPoint, normedDiam, c1, c2));
+                qvartemp.push(Membrane.initial(c1, c2, stats));
+            }
+
+            let jmax = (params.minThickness / params.stepSize) - 1;
+
+            for (let width = 0, widthl = params.maxThickness; width < widthl; width = (jmax + 1) * params.stepSize) {
+                const imax = qvartemp.length - 1 - jmax;
+
+                for (let i = 0, il = imax; i < il; i++) {
+                    const c1 = qvartemp[i].planePoint1;
+                    const c2 = qvartemp[i + jmax].planePoint2;
+
+                    let hphob = 0;
+                    let hphil = 0;
+                    for (let j = 0; j < jmax; j++) {
+                        const ij = qvartemp[i + j];
+                        if (j == 0 || j == jmax - 1) {
+                            hphob += 0.5 * ij.stats.hphob;
+                            hphil += 0.5 * ij.stats.hphil;
+                        } else {
+                            hphob += ij.stats.hphob;
+                            hphil += ij.stats.hphil;
+                        }
+                    }
+
+                    const stats = HphobHphil.of(hphob, hphil);
+
+                    if (hphob > 0) {
+                        const qvaltest = qValue(stats, initialStats);
+                        if (qvaltest > qmax) {
+                            console.log(qmax + ' > ' + qvaltest);
+                            qmax = qvaltest;
+                            membrane = Membrane.scored(c1, c2, HphobHphil.of(hphob, hphil), spherePoint, qmax);
+                        }
+                    }
+                }
+                jmax++;
+            }
+        }
+
+        return membrane!;
+    }
+
+    function qValue(currentStats: HphobHphil, initialStats: HphobHphil): number {
+        if(initialStats.hphob < 1) {
+            initialStats.hphob = 0.1;
+        }
+
+        if(initialStats.hphil < 1) {
+            initialStats.hphil += 1;
+        }
+
+        return (currentStats.hphob * (initialStats.hphil - currentStats.hphil) - currentStats.hphil * (initialStats.hphob - currentStats.hphob)) /
+                (currentStats.total * initialStats.hphob * initialStats.hphil * (initialStats.hphob + initialStats.hphil - currentStats.total));
+    }
+
+    function isInMembranePlane(testPoint: Vec3, normalVector: Vec3, planePoint1: Vec3, planePoint2: Vec3): boolean {
+        const d1 = -Vec3.dot(normalVector, planePoint1);
+        const d2 = -Vec3.dot(normalVector, planePoint2);
+        const d = -Vec3.dot(normalVector, testPoint);
+        return d > Math.min(d1, d2) && d < Math.max(d1, d2);
+    }
+
+    // generates a defined number of points on a sphere with radius = extent around the specified centroid
+    function generateSpherePoints(numberOfSpherePoints: number, centroid: Vec3, extent: number): Vec3[] {
+        const points = new Array<Vec3>(numberOfSpherePoints);
+
+        let oldPhi = 0, h, theta, phi;
+        for(let k = 1, kl = numberOfSpherePoints + 1; k < kl; k++) {
+            h = -1 + 2 * (k - 1) / (numberOfSpherePoints - 1);
+            theta = Math.acos(h);
+            phi = (k == 1 || k == numberOfSpherePoints) ? 0 : (oldPhi + 3.6 / Math.sqrt(numberOfSpherePoints * (1 - h * h))) % (2 * Math.PI);
+
+            const point = Vec3.create(
+                extent * Math.sin(phi) * Math.sin(theta) + centroid[0],
+                extent * Math.cos(theta) + centroid[1],
+                extent * Math.cos(phi) * Math.sin(theta) + centroid[2]
+            );
+            points[k - 1] = point;
+            oldPhi = phi;
+        }
+
+        return points;
+    }
+
+    // generates sphere points close to that of the initial membrane
+    function findProximateAxes(membrane: Membrane, numberOfSpherePoints: number, centroid: Vec3, extent: number): Vec3[] {
+        const points = generateSpherePoints(30000, centroid, extent);
+        const sorted = points.sort((v1, v2) => Vec3.squaredDistance(v1, membrane.point!) - Vec3.squaredDistance(v2, membrane.point!));
+        return sorted.splice(0, numberOfSpherePoints);
     }
 
-    export const enum Flag {
-        NA = 0x0,
-        Membrane = 0x1,
-        NotMembrane = 0x2
+    interface HphobHphil {
+        hphob: number, 
+        hphil: number,
+        total: number
     }
 
-    export function getValue(location: StructureElement.Location, topology: Topology) {
-        const { getSerialIndex } = location.structure.root.serialMapping;
-        const { exposure, serialResidueIndex } = topology;
-        const rSI = serialResidueIndex[getSerialIndex(location.unit, location.element)];
-        if (rSI === -1) return -1;
-        return exposure[rSI];
+    namespace HphobHphil {
+        export function of(hphob: number, hphil: number) {
+            return {
+                hphob: hphob,
+                hphil: hphil,
+                total: hphob + hphil 
+            }
+        }
+
+        const testPoint = Vec3();
+        export function filtered(offsets: ArrayLike<number>, exposed: ArrayLike<boolean>, structure: Structure, label_comp_id: StructureElement.Property<string>, filter: (test: Vec3) => boolean): HphobHphil {
+            const { x, y, z } = StructureProperties.atom;
+            let hphob = 0;
+            let hphil = 0;
+            for (let k = 0, kl = offsets.length; k < kl; k++) {
+                // ignore exposed residues
+                if (exposed[k]) {
+                    continue;
+                }
+    
+                setLocation(l, structure, offsets[k]);
+                Vec3.set(testPoint, x(l), y(l), z(l));
+
+                // testPoints have to be in putative membrane layer
+                if (!filter(testPoint)) {
+                    continue;
+                }
+
+                if (isHydrophobic(label_comp_id(l))) {
+                    hphob++;
+                } else {
+                    hphil++;
+                }
+            }
+            return of(hphob, hphil);
+        }
+    
+        // ANVIL-specific (not general) definition of membrane-favoring amino acids
+        const HYDROPHOBIC_AMINO_ACIDS = ['ALA', 'CYS', 'GLY', 'HIS', 'ILE', 'LEU', 'MET', 'PHE', 'SER', 'THR', 'VAL'];
+        function isHydrophobic(label_comp_id: string): boolean {
+            return HYDROPHOBIC_AMINO_ACIDS.indexOf(label_comp_id) !== -1;
+        }
     }
 
-    export function getFlag(location: StructureElement.Location, topology: Topology) {
-        const value = getValue(location, topology);
-        return value === -1 ? Flag.NA :
-            value ? Flag.NotMembrane :
-                Flag.Membrane;
+    function setLocation(l: StructureElement.Location, structure: Structure, serialIndex: number) {
+        l.structure = structure;
+        l.unit = structure.units[structure.serialMapping.unitIndices[serialIndex]];
+        l.element = structure.serialMapping.elementIndices[serialIndex];
+        return l;
     }
 }

src/tests/browser/render-structure.ts

@@ -148,6 +148,7 @@ async function init() {
         ballAndStick: false,
         molecularSurface: false,
         gaussianSurface: false,
+        membrane: true
     };
 
     const cartoonRepr = getCartoonRepr();
@@ -155,6 +156,7 @@ async function init() {
     const ballAndStickRepr = getBallAndStickRepr();
     const molecularSurfaceRepr = getMolecularSurfaceRepr();
     const gaussianSurfaceRepr = getGaussianSurfaceRepr();
+    const membraneRepr = getMembraneRepr();
 
     if (show.cartoon) {
         cartoonRepr.setTheme({
@@ -201,11 +203,22 @@ async function init() {
         console.timeEnd('gaussian surface');
     }
 
+    if (show.membrane) {
+        // membraneRepr.setTheme({
+        //     color: reprCtx.colorThemeRegistry.create('uniform', { structure }),
+        //     size: reprCtx.sizeThemeRegistry.create('physical', { structure })
+        // });
+        // console.time('membrane layer');
+        // await membraneRepr.createOrUpdate({ ...MembraneRepresentationProvider.defaultValues, quality: 'custom', alpha: 1.0, flatShaded: true, doubleSided: true, resolution: 0.3 }, structure).run();
+        // console.timeEnd('membrane layer');
+    }
+
     if (show.cartoon) canvas3d.add(cartoonRepr);
     if (show.interaction) canvas3d.add(interactionRepr);
     if (show.ballAndStick) canvas3d.add(ballAndStickRepr);
     if (show.molecularSurface) canvas3d.add(molecularSurfaceRepr);
     if (show.gaussianSurface) canvas3d.add(gaussianSurfaceRepr);
+    if (show.membrane) canvas3d.add(membraneRepr);
     canvas3d.requestCameraReset();
     // canvas3d.setProps({ trackball: { ...canvas3d.props.trackball, spin: true } })
 }