In a probit model the data would typically consist of a list of cases, and for each case one would have observed the values of $X_1,X_2,X_3,X_4$ (the number of such predictors may vary) and whether there was a success or a failure. Let $n$ be the number of cases. One observes $(I_i,x_{1,i},\ldots,x_{4,i})$ for $i=1,\ldots,n$, and the value of $I_i$ would be "success" or "failure", and often this would be represented by $1$ or $0$, the number of successes in the $i$th case. But below I'll do it the other way, for sake of consistency with your notation, i.e.
\begin{align}
0 & = \text{success}, \\ 1 & = \text{failure}.
\end{align}
One estimates $\beta_1,\beta_2,\beta_3,\beta_4$ by maximum likelihood. That means the estimates $\widehat\beta_1,\widehat\beta_2,\widehat\beta_3,\widehat\beta_4$ would be the values that maximize the probability of the data that were actually observed. To be more explicit: The likelihood function is
$$
L(\beta_1,\beta_2,\beta_3,\beta_4) = \left(\prod_{i\,:\,I_i=1} \int_{-\infty}^{\beta_1 x_{1,i}+\cdots+\beta_4 x_{4,i}} \varphi(k)\,dk \right)\left(\prod_{i\,:\,I_i=0} \left(1 - \int_{-\infty}^{\beta_1 x_{1,i}+\cdots+\beta_4 x_{4,i}} \varphi(k)\,dk\right) \right).
$$
Each integral is from $-\infty$ to $\beta_1 x_{1,i}+\cdots+\beta_4 x_{4,i}$.
The estimates $\widehat\beta_1,\widehat\beta_2,\widehat\beta_3,\widehat\beta_4$ are the values that maximize this. These would be found by numerical methods; I don't think there's a closed form.
I think typically there would be an additional parameter: one would have $\alpha+\beta_1 x_{1,i}+\cdots+\beta_4 x_{4,i}$ as the upper bound of integration, and the likelihood function would be $L(\alpha,\beta_1,\beta_2,\beta_3,\beta_4)$.
So: No, you don't regress twice.
