I ended up figuring it out and just wanted to post the answer for anyone who may be interested.
Goal
We want to show that linearized gravity predicts that the motion of masses produces gravitational effects similar to those of electromagnetism. Namely, to linear order in the velocity of a test body, the geodesic equation yields an analog of the Lorentz Force equation:
$$\mathbf{a} = - \mathbf{E} - 4\mathbf{v} \times \mathbf{B}.$$
Assumptions
These assumptions are given on page 74 of Wald, but I will summarize them here.
Metric
We are approximating gravity as "weak", so we assume that our metric is given by
$$g_{ab} = \eta_{ab} + \gamma_{ab}, \;\;\;\;\;\;(1)$$
where $\eta_{ab}$ is our flat metric and $\gamma_{ab}$ is a small perturbation of this metric.
We will also be raising and lowering indices with our flat metric $\eta_{ab}$ and $\eta^{ab}$ instead of $g_{ab}$ and $g^{ab}$.
"Linearized" Gravity
For our computations, we will only keep terms that are of linear order in $\gamma_{ab}$, assuming that higher order terms will be sufficiently small.
Covariant Derivative $\nabla_a$
In our first order approximations, we can write $\partial_a$ in place of $\nabla_a$ if applied to $\gamma_{ab}$ or an expression of $\gamma_{ab}$. This is because the Christoffel symbol terms will be not be first order in $\gamma_{ab}$.
Christoffel Symbols
Because $\partial_a \eta_{bc} = 0$, when we write the Christoffel symbols in terms of our flat metric plus our perturbation $\gamma_{ab}$, we get
$$\Gamma^c_{ab} = \frac{1}{2}\eta^{cd}(\partial_a \gamma_{bd} + \partial_b \gamma_{ad} - \partial_d \gamma_{ab}). \;\;\;\;\;\;(2)$$
Gamma Bar $\overline{\gamma}_{ab}$
To simplify computation, we define
$$\overline{\gamma}_{ab} = \gamma_{ab} - \frac{1}{2} \eta_{ab} \gamma, \;\;\;\;\;\;(3)$$
where $\gamma$ is defined as $\gamma_a^a$ as usual.
Note that taking the trace on both sides of (3) results in:
$$\overline{\gamma}_a^{\,\,a} = \gamma_a^{\,\,a} - \frac{1}{2}\eta_a^{\,\,a}\gamma$$
$$\Rightarrow \overline{\gamma} = \gamma - \frac{1}{2} (4) \gamma = -\gamma.$$
Plugging back into (3) gives
$$\gamma_{ab} = \overline{\gamma}_{ab} - \frac{1}{2} \eta_{ab} \overline{\gamma}. \;\;\;\;\;\;(4)$$
Einstein's Equation in the Analog of the Lorentz gauge
Because linearized gravity has a gauge freedom (read page 75 of Wald for more details), we can make a gauge transformation to simplify Einstein's Equation to
$$\partial^c \partial_c \overline{\gamma}_{ab} = -16\pi T_{ab}, \;\;\;\;\;\;(5)$$
for stress-energy tensor $T_{ab}$.
Stress-Energy Tensor
We assume that only lower order effects of motion are taken into account and neglect stresses. Thus, our stress-energy tensor $T_{ab}$ can be approximated to linear order in velocity $t^a$ by
$$T_{ab} = 2t_{(a}J_{b)} - \rho t_a t_b, \;\;\;\;\;\;(6)$$
where $J_b = -T_{ab}t^a$ is the mass current density 4-vector.
Assuming coordinates, for small velocities we can approximate $t^{\mu} = (1,0,0,0)$, so in "matrix" form, our stress-energy tensor can be thought of as
$$T_{\mu \nu} = \begin{bmatrix} \rho & p_x & p_y & p_z \\
p_x & 0 & 0 & 0 \\
p_y & 0 & 0 & 0 \\
p_z & 0 & 0 & 0
\end{bmatrix},$$
where $\mathbf{p} = (p_x, p_y, p_z)$ is momentum.
Analog of Maxwell's Equations
Looking at just the time-time component and space-time components of Einstein's equation, we now have
$$\partial^a \partial_a \overline{\gamma}_{0 \mu} = 16\pi J_{\mu}. \;\;\;\;\;\;(7)$$
Comparing to Maxwell's first equation of E&M in flat 4d-spacetime
$$\partial^a \partial_a A_b = -4\pi J_b,$$
we define
$$A_a \equiv -\frac{1}{4} \overline{\gamma}_{ab} t^b, \;\;\;\;\;\;(8)$$
or for small velocities where $t^{\mu} \approx (1,0,0,0)$,
$$A_{\mu} = -\frac{1}{4} \overline{\gamma}_{0\mu}, \;\;\;\;\;\;(9)$$
We can now define an analog of $\mathbf{E}$ and $\mathbf{B}$ in terms of $A_a = (-\phi, \mathbf{A})$:
$$\mathbf{E} = -\nabla \phi - \frac{\partial \mathbf{A}}{\partial t} \;\;\;\;\;\;(10)$$
$$\mathbf{B} = \nabla \times \mathbf{A}. \;\;\;\;\;\;(11)$$
Derivation
Geodesic Equation
The motion of our test body is given by the geodesic equation,
$$\frac{d^2 x^{\mu}}{d\tau^2} + \Gamma^{\mu}_{\rho \sigma} \frac{d x^{\rho}}{d\tau} \frac{d x^{\sigma}}{d\tau} = 0, \;\;\;\;\;\;(12)$$
where $x^{\mu}(\tau)$ is the world line of the particle in global inertial coordinates. Assuming the motion of our test body is much slower than the speed of light, we approximate proper time $\tau$ by the coordinate time $t$. So, we may write
$$\frac{d^2 x^{\mu}}{dt^2} + \Gamma^{\mu}_{\rho \sigma} \frac{d x^{\rho}}{dt} \frac{d x^{\sigma}}{dt} = 0. \;\;\;\;\;\;(13)$$
Denoting ${d x^{\alpha}}/{dt}$ as $v^{\alpha}$ and moving the second term of (13) to the right side, we get
$$\frac{d^2 x^{\mu}}{dt^2} = -\Gamma^{\mu}_{\rho \sigma} v^{\rho} v^{\sigma}. \;\;\;\;\;\;(14)$$
Lastly, changing $\mu$ to $i$, where $i = 1, 2, 3$, and denoting $a^i = {d^2 x^i}/{dt^2}$, we have
$$a^i = -\Gamma^{i}_{\rho \sigma} v^{\rho} v^{\sigma}, \;\;\;\;\;\;(15)$$
where it is important to note that while $i$ only runs along the spatial indices, $\rho$ and $\sigma$ still run from 0 to 3.
Christoffel Symbols
Approximating time derivatives to be 0, we now calculate our nonzero Christoffel Symbols.
From (2),
$$\Gamma^i_{00} = \frac{1}{2}\eta^{ij}(-\partial_j \gamma_{00})$$
$$=-\frac{1}{2}\partial^i \gamma_{00}.$$
From (4) and (9) and recalling that $A_0 = -\phi$, we find
$$\gamma_{00} = \overline{\gamma}_{00} - \frac{1}{2}\eta_{00 }\overline{\gamma}$$
$$=-4A_0 + \frac{1}{2}(4A_0)$$
$$=4\phi - 2\phi$$
$$=2\phi.$$
So,
$$\Gamma^i_{00} = -\partial^i \phi. \;\;\;\;\;\;(16)$$
Using (2) again and (9), we also have
$$\Gamma^i_{0j} = \Gamma^i_{j0} = \frac{1}{2}\eta^{ik}(\partial_j \gamma_{0k} - \partial_k \gamma_{0j})$$
$$= \frac{1}{2}(\partial_j \gamma_0^i - \partial^i \gamma_{0j})$$
$$= \frac{1}{2}(-4)(\partial_j A^i - \partial^i A_j)$$
$$= 2(\partial^i A_j - \partial_j A^i).$$
Thus, using the fact that $F_{ab} = \partial_a A_b - \partial_b A_a$, we get
$$\Gamma^i_{0j} = \Gamma^i_{j0} = 2F^i_{\,\,j}. \;\;\;\;\;\;(17)$$
Finishing the Computation
Recall that $v^{\alpha} = dx^{\alpha}/{dt}$, so $v^0 = 1$.
Therefore, plugging our results from (16) and (17) into (15) gives
$$a^i = -\Gamma^{i}_{\rho \sigma} v^{\rho} v^{\sigma}$$
$$= -\Gamma^i_{00} - \Gamma^i_{0j}v^j - \Gamma^i_{j0}v^j$$
$$\Rightarrow a^i = \partial^i \phi - 4F^i_{\,\,j}v^j. \;\;\;\;\;\;(18)$$
Writing (10) and (11) in index notation and approximating the time derivative of $\mathbf{A}$ to be 0, we get
$$E^i = -\partial^i \phi \;\;\;\;\;\;(19)$$
$$B^i = \epsilon^{ijk}\partial_j A_k, \;\;\;\;\;\;(20)$$
where $\epsilon^{ijk}$ is the fully antisymmetric Levi Civita symbol (detailed on pages 432 and 433 of Wald).
Lastly, we write $\mathbf{v} \times \mathbf{B}$ in index notation and find
$$(\mathbf{v} \times \mathbf{B})_i = \epsilon_{ijk} v^j B^k$$
$$= \epsilon_{ijk} v^j \epsilon^{klm}\partial_l A_m$$
$$= (\epsilon^{klm} \epsilon_{kij}) v^j \partial_l A_m$$
$$= 2\delta^{[l}_i \delta^{m]}_j v^j \partial_l A_m$$
$$= 2\frac{1}{2}(\delta^{l}_i \delta^{m}_j v^j - \delta^{m}_i \delta^{l}_j v^j) \partial_l A_m$$
$$=v^m \partial_i A_m - v^l \partial_l A_i.$$
So,
$$(\mathbf{v} \times \mathbf{B})_i = F_{ij}v^j. \;\;\;\;\;\;(21)$$
Finally, plugging in our results from (19) and (21) into (18) gives
$$\mathbf{a} = - \mathbf{E} - 4\mathbf{v} \times \mathbf{B}. \;\;\;\;\;\;(22)$$
And we're done.
