Great question. The first part of the answer by niels nielsen is quite inaccurate, talking about "dragging" and "wiggling" of field lines to make the concept more intuitive, but in my opinion it obscures things much more. Zaaikort's answer is great, I only want to add an answer with more detail, and quotes from the reference mentioned in the comments and in Zaaikort's answer.
This explanation is from "Modeling the magnetic pickup of an electric guitar", American Journal of Physics 77, 144–150 (2009). Link.
A simple magnetic pickup (see Fig. 1) is comprised of a permanent magnet surrounded by a coil of wire with typically several thousand turns. The guitar strings consist of wires made of a ferromagnetic material and are parallel to the face of the magnet. The magnetic field of the pickup/wire system, and thus the magnetic flux through the coil, depends critically on the position of the wire. Therefore, moving the wire changes the magnetic flux through the coil. According to Faraday’s law the current induced inside the coil is proportional to the time rate of change of the magnetic flux through the coil. As the string moves through the magnetic field a time-varying current is produced in the coil. This current is used to produce a potential drop across a resistor, which is then amplified and sent to a speaker.
This verifies what you already describe, but then you ask
...how should I model the way the moving steel string is interacting with the magnetic field from the permanent magnet? ... Does it have something to do with the mechanical movement of the electrons in the guitar string as the string vibrates?
First, a tip: Unless you are dealing with quantum effects, you can almost always ignore the role of electrons in electromagnetics. Continuous current distributions are the language of classical E&M, like a material continuum is the language of fluid dynamics. That's at least my opinion, others may disagree.
About the modeling: In short, a ferromagnetic material (the guitar string) has high permeability $\mu$, and thus the resulting field is more intense inside the volume. This is magnetization. The magnetic field around the magnet is not uniform, so changing the position of the string alters the induced magnetic field in the ferromagnetic material, and therefore the total magnetic field is changed (superimpose the field of the magnet and the induced field of the ferromagnetic material).
That's the gist. From an engineering or physics perspective, this situation can be modeled as either (a) a magnetostatic problem, (b) a magnetic circuit, or (c) a problem of (fictitious) magnetic charge density. Each of the representations are suitable, they represent the same phenomena, just at different levels of generality and "physicality". To develop a model, first begin with a description of the magnetic field due to the permanent magnet, without the guitar string. Most models of permanent magnets are too complicated to be useful outside of computation, but using the fictitious concept of magnetic charge yields a good approximation with simpler expressions, and this is the route taken by the above-cited paper. Then, determine the effect of the nearby ferromagnetic material. This might involve assuming the wire to be a series of parallel plate magnets with no width and a height equal to the diameter of the wire, above the permanent magnet. From the paper:
... We model the wire as a series of infinitesimally wide magnets whose strength is linearly proportional to the local field at the position of the wire due to the permanent magnet, and the height of these infinitesimal magnets is equal to the diameter of the wire.
We assume that the coercitive field of the permanent magnet is sufficiently large so that the presence of the wire does not affect the magnetic charge density on the surface. Therefore, the only change in the magnetic field at $(x_p,y_p,z_p)$ due to the presence of the wire is caused by the magnetization of the wire itself. We further assume that the magnetic intensity at the wire is less than is required for saturation, and we ignore the effects of hysteresis. Therefore, the magnetic field of the permanent magnet results in a linear change in magnetization at the wire.
Finally, the model itself is calculated from what is essentially Coulomb's law, but for magnetic charge. It is calculated in three parts. First, the field due to the permanent magnet; then, the magnitude of the magnetic field at a given position is calculated; finally, the z-component of the magnetic field is used to calculate flux.
I will reproduce the reasoning in the linked paper for completeness. Let the radius of the permanent magnet be $\psi$, with center $(x_0,y_0,z_0)$, with its face perpendicular to the $z$-axis, and plate separation equal to its length. Denote the magnetic charge density by $\sigma$, and assume this density is uniform. Then the $z$-component of the magnetic field at $(x',y',z')$ is
$$
B_{z}(x',y',z') = \int_{0}^{2\pi} \int_{0}^{\psi} \frac{\sigma(z'-z_0)\rho}{\left[ (x'-[x_0-\rho \cos \phi])^{2} + (y' - [y_0-\rho \sin \phi]) ^{2} + (z'-z_0)^{2} \right]^{3/2} } d\rho d\phi,
$$
Next, modeling the wire as an infinitesimally wide disk located at $(x',y',z')$, the magnitude of the magnetic field at this disk, before taking the permeability into account, is:
$$
|\mathbf {B}_w (x',y',z')| = \int_{0}^{2\pi} \int_{0}^{\psi} \frac{\sigma\rho}{\left[ (x'-[x_0-\rho \cos \phi])^{2} + (y' - [y_0-\rho \sin \phi]) ^{2} + (z'-z_0)^{2} \right] } d\rho d\phi,
$$
(Note: subscript $w$ is for "wire") Finally, the z-component of this field is calculated by multiplying the above expression with a constant of proportionality related to the permeability/susceptibility of the wire, and resolving it into its z-component, for an arbitrary position $(x,y,z)$, which is the location of interest for the total magnetic field,
$$
B_{w,z}(x,y,z) = \gamma |\mathbf B_w(x',y',z')| \frac{ (z'-z) }{\left[(x'-x)^2 + (y'-y)^2 +(z'-z)^2\right]^{3/2}},
$$
and this is summed for each segment of the wire. In other words, splitting the wire into $N$ pieces, with positions $(x_i,y_i,z_i)$ for $i=1,2,\cdots,N$, the field at $(x,y,z)$ is
$$
B_{w,z}^{tot}(x,y,z) = \gamma \sum_i |\mathbf B_w(x_i,y_i,z_i)| \frac{ (z_i-z) }{\left[(x_i-x)^2 + (y_i-y)^2 +(z_i-z)^2\right]^{3/2}}.
$$
