Modifying Spline Algorithms to Draw Proteins with WebAssembly

This post is a modified version of a presentation I gave about my protein visualiser available online with WebAssembly

Proteins are the building blocks of our cells and allow for the incredibly complicated behaviours of the body. Knowing how important proteins are then, models of protein’s structures allow scientists to quickly compare and understand a protein’s features. One protein model is the ribbon diagram which has been around since the 1980s. They were originally hand drawn but were quickly created by computers. A Ribbon Diagram shows the shape of a protein by drawing ‘ribbons’ between the alpha-carbons (which are central to the peptide bonds of a polypeptide chain) of each amino-acid. Looking at only the alpha-carbons simplifies the thousands of individual atoms down to relatively few points which can be connected into a shape.

Before we can create our diagrams we first need the atoms which make up our protein, thankfully there is a well established (and documented!) format (PDB) along with a database of proteins called the Worldwide Protein Data Bank. While the format is a bit old, it really isn’t too bad, my reader (See loadPdbFile) came out at ~300 lines of . Or you could use one of many libraries suggested by the PDB people. Either way you now have your atoms and bonds loaded into a nice and pretty data structure. For my online viewer, reading is quick enough that I don’t even bother caching anything, the PDB is processed every page load although this could certainly be improved. On to rendering!

To draw our proteins, like pretty much all 3D graphics, we need to form meshes which the GPU can understand. These meshes consist of triangles, each consisting of three vertices in 3D space which together enclose the ribbon of a protein, for correct lighting we also need a normal for each vertex. I’m going to skip over the details of my renderer, ultimately you could draw the meshes with any 3D program you liked or even export 3D models.

Creating a ribbon consists of two steps, fitting a spline to the alpha-carbons and extruding a ‘profile’ along this spline. A spline takes a pair of points and creates a curve which consists of many points which smoothly move between the pair. To extrude our spline we have to know the normal and tangent at each of these points. A profile is the shape of the ribbon’s cross section and can be circular or rectangular.

Before creating the spline I first average adjacent alpha-carbon positions to create a smoother look, we then need to choose tangents (parallel to tube) for each point. You could use more chemically accurate tangents but I used a Cardinal interpolation to achieved a smoother look. The Cardinal method uses a dampened vector between adjacent position, so that the tangent is given by

Now we need to actually generate the spline, it’s important that our spline not only smoothly interpolates position but also tangent, normals and binormals. That is, the spline’s first and second derivatives must be continuous. Specifying the tangents at each point ensures that adjacent splines have matching first derivatives but not second derivates. If we want we could use a quadratic spline, but we sacrifice smoothness and I found that this produces a poor result. Rather, I discard the true normal and calculate the normal from the current tangent and previous point’s normal. For smoothness set the endpoint second derivative to zero, as in the natural Hermite Cubic. This ensures continuity, since the continuity of the tangents are guaranteed and simplifies calculation. If the previous normal is invalid I set the binormal to the perpendicular component, relative to the tangent, of the average of the endpoints’ binormals, but I’m sure there are better choices.

The position , where , is given by

And the tangent is

The tangent is then be normalized and used to calculate the binormal and normal

Or in code if you prefer (simplified to remove profile blending and doubled normals)

auto bref = (b0 + b1) / 2.f;
for(int i = 0; i < num_points_per_spline; i++) {
    auto t1 = (float)i / (num_points_per_spline-1);
    auto t2 = t1*t1;
    auto t3 = t1*t2;

    auto p = (2*t3-3*t2+1)*pp0 
           + (t3-2*t2+t1)*m0
           + (-2*t3+3*t2)*pp1
           + (t3-t2)*m1;

    auto tn = (6*t2-6*t1)*pp0 
            + (3*t2-4*t1+1)*m0
            + (-6*t2+6*t1)*pp1
            + (3*t2-2*t1)*m1;

    tn = normalize(tn);
    vec3 bn;
    if(length(prev_normal) > EPSILON) {
        bn = normalize(cross(tn, prev_normal));
    } else {
        bn = normalize(perpendicularComponent(bref, tn));
    }
    auto n = cross(bn, tn);

    projectPointsOnPlane(num_splines, p, bn, n, pf, &points[num_splines*i]);
    projectPointsOnPlane(num_splines, vec3(0), bn, n, pfn, &normals[num_splines*i]);
    prev_normal = n;
}