We are going to create a single-page web app that runs a Gray-Scott reaction-diffusion simulation.
 The simulation will run on a canvas and display organic Turing patterns in real time.
 The mouse will be used to seed new growth by adding chemicals at the mouse position.

 Requirements:
 - All code in one file: index.html
 - CSS in a <style> tag, JavaScript in a <script> tag
 - No external resources, everything inline

 Approach:
 1. We'll create an HTML5 canvas element.
 2. We'll use WebGL for the simulation because it's computationally intensive and WebGL (via fragment shaders) can run efficiently on the GPU.
 3. The Gray-Scott model involves two chemicals (U and V) that react and diffuse. We'll simulate using a fragment shader that does the computation.
 4. We'll use a ping-pong rendering technique: two textures (or framebuffers) that we swap each frame to use the previous state as input for the next.
 5. The mouse interaction will add a spot of V (the second chemical) at the mouse position.

 Steps in the code:
 - Set up the WebGL context and shaders.
 - Create two textures and framebuffers for the ping-pong rendering.
 - Write the fragment shader that computes the Gray-Scott reaction-diffusion.
 - Write the vertex shader and a simple fragment shader for rendering the result to the screen.
 - Handle mouse events to add V at the mouse position.

 Gray-Scott parameters (example values that produce interesting patterns):
   f: feed rate (e.g., 0.055)
   k: kill rate (e.g., 0.062)
   dU: diffusion rate for U (e.g., 0.16)
   dV: diffusion rate for V (e.g., 0.08)

 However, note: we'll set the diffusion rates in the shader and allow interaction to change the pattern.

 Implementation:

 We'll have:
   - A vertex shader for a full-screen quad (simple).
   - A fragment shader for the simulation step (Gray-Scott update).
   - A fragment shader for displaying the result (which might just take the V value and color it).

 We'll use two programs:
   Program 1: simulation program (uses the Gray-Scott shader) to update the state.
   Program 2: rendering program to draw the current state to the screen.

 We'll use two framebuffers (with attached textures) to hold the state (U and V for each pixel). We can store U and V in two channels of a texture (e.g., R and G).

 Steps each frame:
   1. Bind framebuffer A (as the render target) and use the simulation program with the state from framebuffer B (as the input texture).
   2. Render a full-screen quad to run the simulation shader on every pixel.
   3. Swap: now framebuffer A has the updated state. We then use framebuffer A as input for the next step.
   4. For rendering to the screen: use the rendering program with the updated state (from framebuffer A) and draw to the canvas.

 Mouse interaction:
   When the mouse is moved (or dragged), we'll update a texture that represents the user input (a spot of V). However, we can also handle it by updating the state texture directly in the simulation step by checking the mouse position.

 Alternatively, we can pass the mouse position and a flag to the simulation shader and have it add V at that position.

 However, note: the simulation runs on the entire texture. We can't write to arbitrary locations without breaking the simulation. Instead, we can do:

   - In the simulation shader, check if the current fragment is near the mouse and then add some V.

 But note: the simulation shader runs for every pixel. We can pass the mouse position and then in the shader, if the distance is small, add to V.

 We'll pass:
   u_mouse: the mouse position in texture coordinates (0 to 1) and whether the mouse is down.

 We'll have to update the state texture in the simulation step accordingly.

 Let's design:

   Uniforms for the simulation shader:
     u_texture: the state texture from the previous frame (with two channels: R=U, G=V)
     u_resolution: the resolution of the canvas
     u_mouse: vec3 (x, y, button) where button is 1 if pressed, 0 otherwise
     u_time: for any time-dependent effects? (not necessary for Gray-Scott, but maybe for changing parameters)

   In the simulation shader, we'll compute the next U and V.

   The Gray-Scott equations:
     u_next = u + (dU * laplacian(u) - u*v*v + f*(1 - u)) * dt
     v_next = v + (dV * laplacian(v) + u*v*v - (f+k)*v) * dt

   We'll compute the laplacian using a 3x3 convolution kernel (e.g., using the Moore neighborhood).

   We'll use a time step (dt) of 1.0, but we might adjust.

   For the laplacian, we can use a kernel like:

     [ 0.05, 0.2, 0.05 ]
     [ 0.2, -1.0, 0.2 ]
     [ 0.05, 0.2, 0.05 ]

   But note: the actual diffusion kernel for Gray-Scott is often a 3x3 with center -1 and neighbors 0.2? We'll use a common one.

   Alternatively, we can use a 9-point stencil:

     laplacian = (neighbor_sum - center * 8) / 1.0   [but then we multiply by diffusion rate]

   Actually, the common kernel for Gray-Scott is:

     kernel = [ [1,2,1], [2,-12,2], [1,2,1] ] / 6.0   ? 

   I found an example that uses:

     laplacian = ( 
        texture2D(u_texture, coord + offset[0]) +
        texture2D(u_texture, coord + offset[1]) +
        ... for all 8 neighbors
        - 8.0 * center ) / 3.0;

   But let's stick with a common one: 

     We'll use the 5-point stencil for simplicity? Actually, let's use 9-point:

     kernel = [0.05, 0.2, 0.05,
               0.2, -1.0, 0.2,
               0.05, 0.2, 0.05]

   Then the laplacian is the sum of the neighbors multiplied by their kernel weights.

   We'll precompute the offsets for the 8 neighbors.

   However, note: we are storing U and V in two channels. We'll compute the laplacian for U and V separately.

   Steps in the simulation shader:

     vec2 texCoord = gl_FragCoord.xy / u_resolution;
     vec2 center = texture2D(u_texture, texCoord).rg;

     // Compute laplacian for U and V
     vec2 laplacian = vec2(0.0);
     for (int i = -1; i <= 1; i++) {
        for (int j = -1; j <= 1; j++) {
           if (i == 0 && j == 0) continue;
           vec2 offset = vec2(i, j) / u_resolution;
           vec2 neighbor = texture2D(u_texture, texCoord + offset).rg;
           // Use a kernel: for simplicity, we'll use the same weight for all neighbors? 
           // Actually, we'll use a weight of 0.2 for the direct neighbors (up, down, left, right) and 0.05 for diagonals.
           float weight = (abs(i) + abs(j) == 1) ? 0.2 : 0.05;
           laplacian += neighbor * weight;
        }
     }
     laplacian -= center * (0.2*4 + 0.05*4);   // because we added 4*0.2 and 4*0.05 to the neighbors, then subtract center*(sum of weights for neighbors) but note: the center weight is negative.

     Actually, the kernel we are using is:

        [ 0.05, 0.2, 0.05 ]
        [ 0.2, -1.0, 0.2 ]
        [ 0.05, 0.2, 0.05 ]

     So the center weight is -1.0? Then:

        laplacian = (neighbor_sum * 1.0) - center * 1.0? 

     But note: the kernel is applied as:

        result = center * (-1.0) + neighbor_sum * (with weights summing to 1.0? Actually, the weights for neighbors sum to 1.0? Then the entire kernel sums to 0).

     Actually, we can compute:

        laplacian = -center;
        for each neighbor, add neighbor * weight.

     Then we have:

        laplacian = (neighbor_sum * weight) - center * (total_weight_of_neighbors) ?

     Alternatively, we can precompute the kernel weights and then:

        laplacian = 0.0;
        for each offset, laplacian += weight * neighbor;
        then subtract center * (sum of weights) ?

     Actually, the kernel is:

        [ w00, w01, w02 ]
        [ w10, w11, w12 ]
        [ w20, w21, w22 ]

        and w11 = -1.0, and the others are positive and sum to 1.0? Then:

        laplacian = (w00*tex00 + w01*tex01 + ... + w22*tex22) 

        and then we use that.

     We'll use:

        float kernel[9] = float[](
            0.05, 0.2, 0.05,
            0.2, -1.0, 0.2,
            0.05, 0.2, 0.05
        );

        Then we traverse the 3x3 grid and accumulate.

     But note: the kernel is applied to the entire 3x3, including the center.

     So:

        vec2 laplacian = vec2(0.0);
        int index = 0;
        for (int i = -1; i <= 1; i++) {
            for (int j = -1; j <= 1; j++) {
                vec2 neighbor = texture2D(u_texture, texCoord + vec2(i,j)/u_resolution).rg;
                float weight = kernel[index];
                index++;
                laplacian += neighbor * weight;
            }
        }

     Then we have the laplacian without subtracting the center separately.

   Then we compute:

        float u = center.r;
        float v = center.g;
        float du = dU * laplacian.r;
        float dv = dV * laplacian.g;

        float reaction = u * v * v;
        float u_new = u + (du - reaction + f * (1.0 - u)) * dt;
        float v_new = v + (dv + reaction - (f + k) * v) * dt;

   Then we clamp to [0,1].

   Additionally, if the mouse is near, we add V:

        if (u_mouse.z > 0.5) { // button pressed
            float dist = distance(gl_FragCoord.xy, u_mouse.xy);
            if (dist < 10.0) { // within 10 pixels
                v_new = min(1.0, v_new + 0.5);
            }
        }

   Then we output: vec4(u_new, v_new, 0.0, 1.0);

   For the rendering shader, we can simply output a color based on V (or U and V). A common way is to use V for color:

        gl_FragColor = vec4(vec3(v), 1.0);

   Or we can use a colormap.

   We'll start with grayscale: V as the color.

 Implementation steps:

   We'll write the HTML structure with a canvas.

   We'll set up WebGL, create two framebuffers with textures.

   We'll compile the shaders:

     - Vertex shader (for both programs): simple pass-through.

     - Simulation fragment shader: does the Gray-Scott update.

     - Render fragment shader: displays the V value.

   We'll create two programs: simulationProgram and renderProgram.

   We'll create a full-screen quad.

   Then in the animation loop:

     - Use simulationProgram to update the state to framebuffer A (using texture B as input).
     - Then swap: use framebuffer A as input for the next frame.
     - Then use renderProgram to draw the state (from framebuffer A) to the screen.

   We'll handle mouse events to update a uniform for the mouse.

   We'll set the parameters f, k, dU, dV as uniforms. We can adjust them to see different patterns.

   We'll initialize the state with a central square of V and the rest U=1, V=0.

   Alternatively, we can initialize with random noise? Or a central spot.

   We'll start with a central square:

        for each pixel:
          if (x,y) is near center, set U=0.5, V=0.25? Or according to Gray-Scott: U=1, V=0 everywhere except a square in the center where V=1.

   Actually, we can initialize in the first frame by drawing to the state texture.

   We'll do:

        Create an initial state texture: U=1.0 everywhere, V=0.0 except in a rectangle in the center: V=1.0.

   But note: we are using two framebuffers. We can initialize both with the same initial state.

   Steps:

        Create textureA and textureB, and framebufferA and framebufferB.

        Initialize textureA: create a Float32Array and set the data for the initial state.

        Then we set textureA as the initial state and textureB as a copy? Actually, we can initialize both the same.

        Then we start the loop: first we render to framebufferB using textureA as input, then swap.

   We'll set the initial state in JavaScript by creating a texture with initial data.

   Alternatively, we can use a shader to initialize? But for simplicity, we'll initialize with JavaScript.

   We'll create a Float32Array for the initial state (for one texture) and then upload to both textures? Actually, we only need one initial state. Then we start by rendering from textureA to framebufferB? Then we swap.

   Actually, we'll start with:

        state1 = initial state
        state2 = ... (we don't care, because we'll update immediately)

        Then:

        frame 0: use state1 to update to state2 (render to framebufferB using state1 as input)
        frame 1: use state2 to update to state1 (render to framebufferA using state2 as input)
        frame 2: use state1 to update to state2, etc.

   So we start by having state1 set to the initial state. Then we render the simulation to framebufferB (state2) using state1. Then we use state2 as input for the next step to update state1.

   Therefore, we only need to initialize state1.

   Steps:

        Create two textures: stateTexture1 and stateTexture2, and two framebuffers: fb1 and fb2.

        Attach stateTexture1 to fb1, stateTexture2 to fb2.

        Initialize stateTexture1 with initial data.

        Then in the loop:

            step 0: bind fb2, use stateTexture1 as input -> output to stateTexture2
            step 1: bind fb1, use stateTexture2 as input -> output to stateTexture1
            step 2: bind fb2, use stateTexture1 as input -> output to stateTexture2
            ...

        Then to render to the screen, we use the stateTexture that was just updated? Actually, we are alternating. We can use the one that was updated in the previous step.

        But note: after each simulation step, we swap the framebuffers and the textures. Then we render the latest state to the screen.

        We'll keep track: currentStateTexture and nextStateFramebuffer.

        Actually, we'll maintain:

            readFramebuffer: the one we read from (the state from the previous frame)
            writeFramebuffer: the one we write to (the next state)

        Then after writing, we swap: 
            temp = readFramebuffer; readFramebuffer = writeFramebuffer; writeFramebuffer = temp;

        But note: the framebuffers are bound to textures. We'll also swap the textures? Actually, we can swap the framebuffers and then the read texture is the one attached to the readFramebuffer.

        Alternatively, we can have two textures and swap them:

            let currentTexture = stateTexture1;
            let nextTexture = stateTexture2;

            Then in the simulation step:
                bind framebuffer that has nextTexture attached.
                set currentTexture as the input texture.
                render.

            Then swap: 
                [currentTexture, nextTexture] = [nextTexture, currentTexture];

        Then for rendering to the screen, we use currentTexture (which is the one we just wrote to? Actually, after the simulation step, we swapped so currentTexture becomes the next state). 

        Actually, after the simulation step:

            currentTexture (was the input) becomes the old state, and nextTexture becomes the new state.

            Then we swap: currentTexture = nextTexture, and nextTexture = currentTexture? 

            Actually, we want to use the new state for rendering? Then after the simulation step, we have the new state in nextTexture. Then we swap so that nextTexture becomes the currentTexture for the next simulation step? And we use the nextTexture (which is now the new state) for rendering? 

            Let me clarify:

                Step 0:
                  currentTexture = stateTexture1 (initial state)
                  nextTexture = stateTexture2 (undefined)

                  We bind framebuffer for nextTexture (fb2) and set currentTexture (stateTexture1) as input -> output to nextTexture (stateTexture2).

                  Then we swap: 
                      let temp = currentTexture; currentTexture = nextTexture; nextTexture = temp;

                  Now currentTexture = stateTexture2 (the new state) and nextTexture = stateTexture1.

                Step 1:
                  Bind framebuffer for nextTexture (fb1) and set currentTexture (stateTexture2) as input -> output to nextTexture (stateTexture1).

                  Then swap: currentTexture = stateTexture1, nextTexture = stateTexture2.

            Then for rendering to the screen, we use currentTexture (which is the state we just computed) right after the simulation step? Actually, we do the simulation step and then we render the currentTexture to the screen? 

            But note: after the simulation step, we have the new state in nextTexture (before swap) and then we swap so that currentTexture becomes the new state. Then we render currentTexture.

        So:

            for each frame:
                bind write framebuffer (which is the framebuffer attached to nextTexture)
                set the simulation program and set the uniform for the state texture to currentTexture.
                render the full-screen quad -> writes to nextTexture.

                swap: [currentTexture, nextTexture] = [nextTexture, currentTexture];

                Then bind the default framebuffer (canvas) and use the render program with currentTexture (which is the one we just wrote to) as input.

                Then request next frame.

   We'll set up:

        let currentTexture = stateTexture1;
        let nextTexture = stateTexture2;

        let currentFramebuffer = fb1;   // not needed, because we bind the framebuffer that has nextTexture attached? Actually, we have two framebuffers: fb1 is attached to stateTexture1, fb2 to stateTexture2.

        So when we want to write to nextTexture, we bind the framebuffer that has nextTexture attached.

        How we know which framebuffer? We can store:

            framebuffers: [fb1, fb2] and textures: [tex1, tex2]

        Then when we write to nextTexture (which is tex2 initially), we bind fb2.

        Then after swap, nextTexture becomes tex1, so we bind fb1.

        Alternatively, we can have an array and swap:

            let readIndex = 0;
            let writeIndex = 1;

            Then:
                bind framebuffer[writeIndex]
                set texture uniform to texture[readIndex]
                render simulation

                then swap: 
                    let temp = readIndex;
                    readIndex = writeIndex;
                    writeIndex = temp;

            Then for rendering to the screen, we use texture[readIndex] (because after swap, readIndex points to the one we just wrote).

   We'll do the array index method.

   Steps:

        let textures = [tex1, tex2];
        let framebuffers = [fb1, fb2];
        let readIndex = 0;
        let writeIndex = 1;

        Then in the draw loop:

            // Bind the framebuffer we are writing to (writeIndex)
            gl.bindFramebuffer(gl.FRAMEBUFFER, framebuffers[writeIndex]);
            gl.viewport(0,0,width,height);
            gl.useProgram(simulationProgram);
            // Set the state texture: textures[readIndex]
            gl.activeTexture(gl.TEXTURE0);
            gl.bindTexture(gl.TEXTURE_2D, textures[readIndex]);
            gl.uniform1i(gl.getUniformLocation(simulationProgram, "u_texture"), 0);

            // Set other uniforms: resolution, mouse, time, etc.

            // Draw the quad

            // Now swap:
            let temp = readIndex;
            readIndex = writeIndex;
            writeIndex = temp;

            // Now bind the canvas framebuffer (null) to draw to screen
            gl.bindFramebuffer(gl.FRAMEBUFFER, null);
            gl.viewport(0,0,width,height);
            gl.useProgram(renderProgram);
            // Set the texture to the one we just updated: textures[readIndex] (because we swapped, readIndex is now the one we wrote to)
            gl.activeTexture(gl.TEXTURE0);
            gl.bindTexture(gl.TEXTURE_2D, textures[readIndex]);
            gl.uniform1i(gl.getUniformLocation(renderProgram, "u_state"), 0);

            // Draw the quad

            requestAnimationFrame(draw);

   We'll initialize the state of texture[0] (tex1) to the initial state.

   Initial state: 

        We create an array of data: for each pixel, 
          U = 1.0
          V = 0.0

        Then in a central rectangle, set V to 1.0.

        We'll use:

          let data = new Float32Array(width * height * 4); // RGBA, we only need RG for U and V, but we have 4 channels
          for (let i=0; i<width*height; i++) {
            let x = i % width;
            let y = Math.floor(i / width);
            data[i*4] = 1.0;   // U
            data[i*4+1] = 0.0; // V
            // set center: 
            if (Math.abs(x - width/2) < 10 && Math.abs(y - height/2) < 10) {
                data[i*4+1] = 1.0;
            }
          }

        Then upload to texture0.

   We'll set the texture parameters for floating point and linear filtering? Actually, we need nearest neighbor for the simulation? Because we are doing discrete convolution.

        gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST);
        gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST);

   We'll set the internal format to RGBA32F and format to RGBA, type to FLOAT.

   We'll set the canvas size to, say, 512x512.

   We'll handle window resize? We'll fix the size for simplicity.

   We'll set the vertex shader for both programs:

        attribute vec2 a_position;
        varying vec2 v_texCoord;
        void main() {
            gl_Position = vec4(a_position, 0.0, 1.0);
            v_texCoord = (a_position + 1.0) / 2.0; // because position is in [-1,1] and we want [0,1] for texCoord
        }

   Then the fragment shaders get v_texCoord.

   For the simulation fragment shader, we use:

        uniform sampler2D u_texture;
        uniform vec2 u_resolution;
        uniform vec3 u_mouse; // (x, y, button)

        varying vec2 v_texCoord;

        void main() {
            // We are working in pixel coordinates: gl_FragCoord.xy
            vec2 position = gl_FragCoord.xy;

            // Parameters (we can make them uniforms too)
            float dU = 0.16;
            float dV = 0.08;
            float f = 0.055;
            float k = 0.062;
            float dt = 1.0;

            // Current state at this pixel
            vec2 center = texture2D(u_texture, v_texCoord).rg;

            // Compute laplacian using 3x3 kernel
            vec2 laplacian = vec2(0.0);
            float kernel[9] = float[](
                0.05, 0.2, 0.05,
                0.2, -1.0, 0.2,
                0.05, 0.2, 0.05
            );
            int index = 0;
            for (int i = -1; i <= 1; i++) {
                for (int j = -1; j <= 1; j++) {
                    vec2 offset = vec2(i, j) / u_resolution;
                    vec2 neighbor = texture2D(u_texture, v_texCoord + offset).rg;
                    laplacian += neighbor * kernel[index];
                    index++;
                }
            }

            // Compute reaction
            float u = center.r;
            float v = center.g;
            float reaction = u * v * v;

            float du = dU * laplacian.r;
            float dv = dV * laplacian.g;

            float u_new = u + (du - reaction + f * (1.0 - u)) * dt;
            float v_new = v + (dv + reaction - (f + k) * v) * dt;

            // Clamp
            u_new = clamp(u_new, 0.0, 1.0);
            v_new = clamp(v_new, 0.0, 1.0);

            // Mouse interaction: if mouse button is down and near the mouse position, add V
            if (u_mouse.z > 0.5) {
                float dist = distance(position, u_mouse.xy);
                if (dist < 10.0) {
                    v_new = min(1.0, v_new + 0.5);
                }
            }

            gl_FragColor = vec4(u_new, v_new, 0.0, 1.0);
        }

   For the render fragment shader:

        uniform sampler2D u_state;
        varying vec2 v_texCoord;

        void main() {
            float v_val = texture2D(u_state, v_texCoord).g;
            gl_FragColor = vec4(vec3(v_val), 1.0);
        }

   We'll also pass the resolution to the simulation shader.

   We'll set up the mouse events to update a uniform in the simulation program.

   We'll create a uniform for u_mouse (vec3) and update it on mousemove and mousedown/up.

   We'll set:

        canvas.addEventListener('mousemove', function(e) {
            var rect = canvas.getBoundingClientRect();
            mouseX = e.clientX - rect.left;
            mouseY = rect.height - (e.clientY - rect.top); // because in WebGL, bottom is 0? Actually, our coordinate: gl_FragCoord.y is from bottom? No, it's from the bottom of the canvas? 

            Actually, gl_FragCoord.y is from the bottom? But in the mouse event, clientY is from the top. So:

            mouseY = rect.height - (e.clientY - rect.top); // to flip to bottom-up?

            But note: our texture coordinates are from bottom? Actually, in the vertex shader we set v_texCoord = (a_position+1)/2, and a_position goes from -1 to 1 with y up? Then the texture is sampled with (0,0) at bottom left? 

            However, when we draw the canvas, the canvas's coordinate system in the browser has (0,0) at top left. But our WebGL drawing uses bottom-left for the texture? 

            Actually, we are rendering to the canvas without flipping. So the canvas has (0,0) at top left? Then we have to flip the mouseY? 

            Alternatively, we can set the canvas to use top-left for the fragment shader? 

            We can set the projection matrix? But we are using a full-screen quad with coordinates from -1 to 1, and the canvas has (0,0) at top left? Then gl_FragCoord.y is from the top? 

            Actually, gl_FragCoord.y: the window space origin is at the lower-left corner. So if we have a canvas of height H, then gl_FragCoord.y=0 at the bottom and H at the top? 

            But the mouse event gives us y from top (0) to bottom (H) at the top? 

            So we need:

                mouseY = H - (e.clientY - rect.top); 

            But note: rect.top is the top of the canvas relative to the viewport. 

            Alternatively, we can set:

                mouseY = rect.height - (e.clientY - rect.top);

            This gives the bottom-up coordinate.

            However, let me test: 

                When the mouse is at the top: e.clientY - rect.top = 0 -> mouseY = rect.height -> which is the top in gl_FragCoord? 

                But gl_FragCoord.y=0 at the bottom and rect.height at the top? 

            So yes, we need to flip.

        });

        And on mousedown, set mouseButton=1, on mouseup set mouseButton=0.

        Then set the uniform: u_mouse = [mouseX, mouseY, mouseButton]

   We'll update this uniform every frame? Or every mousemove? We can update it in the draw loop? Because the mouse might move without the frame updating? But we update every frame, so we can update the uniform in the draw loop.

   We'll store the current mouse state in variables and update the uniform in the simulation program at the beginning of the draw loop.

   We'll set the uniform location for u_mouse and then set it.

   We'll write the code accordingly.

   Let's code:

        We'll create the HTML file with a canvas and the script.

   Note: We are writing everything in one file.

   We'll set the canvas size to 512x512.

   We'll include the WebGL context creation.

   We'll handle WebGL errors.

   We'll write the shaders as strings.

   We'll compile and link the programs.

   We'll create the full-screen quad.

   We'll initialize the state texture.

   We'll set up the animation loop.

   We'll handle mouse events.

   We'll write the shader code as multiline strings.

   We'll use highp precision for the fragment shaders.

   Let's write the HTML structure:

        <!DOCTYPE html>
        <html>
        <head>
            <title>Gray-Scott Reaction-Diffusion</title>
            <style>
                body { margin: 0; }
                canvas { display: block; }
            </style>
        </head>
        <body>
            <canvas id="canvas" width="512" height="512"></canvas>
            <script>
                // All JavaScript code here
            </script>
        </body>
        </html>

   We'll write the JavaScript code inside the script tag.

   Due to the complexity, we'll break the code into sections.

   Note: We are not using any external libraries.

   Let's start coding.

   We'll write the entire index.html file.

   We'll write the shader strings.

   We'll create the WebGL context.

   We'll create the two textures and framebuffers.

   We'll initialize the state texture.

   We'll set up the programs.

   We'll set up the quad.

   We'll set up the mouse handling.

   We'll write the draw loop.

   We'll set the parameters.

   We'll start the animation.

   This code will be long, so we'll write it step by step.

   Note: We are writing to index.html in the current directory.

   We'll write the entire content in one go.

   Let's write the file.